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SPARSE  MATRIX  SYMPOSIUM 

Fairfield  Glade,  Tennessee 
October  24-27,  1902 


Sunday,  October  24 

6:00-10:00  p.m.  Registration. 

8:00-10:00  p.m.  Reception. 


Monday,  October  25** 

8:45  a.m.  Welcoming  Remarks,  Robert  C.  Ward. 

Session  1.  Invited  Papers.  St.  George  Room 
Chairperson:  Robert  C.  Ward 

9:00  a.m.  Iain  S.  Duff*  and  John  K.  Reid,  "Direct  Methods  for 
Solving  Sparse  Systems  of  Linear  Equations." 

9:45  Stanley  C.  Eisenstat,  "Iterative  Methods  for  Solving  Large 

Sparse  Linear  Systems." 


10:30  Coffee  break. 

Session  2A.  Contributed  Papers.  Dru id-Dorchester  Room 
Chairperson:  J.  4 lan  George 

11:00  a.m.  Albert  M.  Erisman,  Roger  G.  '..rimes,  John  G.  Lewis*,  a  no 
William  G.  Poole,  !r. ,  A  at-  u. rurally  stable  Modifi¬ 
cation  of  the  Mel  1  erma.n-Rari  Algorithm  for  Reorder  m 
Unsymmetric  Sparse  Matrices. 

11:20  Arne  Drud,  "Spi  vp  select  ion  f  .  arge  f»r.arse  Sets  o 

Nonlinear  Equations. 

11:40  Robert  Fourer,  "  So  *  v  i  n Sparse  stair-,  ase  Systems  by 

Gauss i an  El imi nat i on . 

♦Presenter  of  paper. 

♦♦Poster  Session  open  during  all  <  1 
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Session 

11:00  a 

11:20 

11:40 

12:00 

Session 

12:10  p 

12:30 

12:50 

Session 

12:10  p 

12:30 

12:50 


Monday,  October  25** 


2B.  Contributed  Papers.  Dru i d-Dorchester  Room 
Chairperson:  Carlos  A.  Felippa 

m.  R.  L.  Cox,  "LS0D28:  A  Variant  of  LS0DE  for  Problems 
Having  a  General  Large  Sparse  Jacobian." 

Granville  Sewell,  "IMSL  Software  for  Differential 
Equations  in  One  Space  Variable." 

J.  W.  Neuberger  and  R.  J.  Renka*,  "A  Portable  Software 
Package  for  Nonlinear  Partial  Differential  Equations." 


Break . 


3A.  Contributed  Papers.  St.  George  Room 
Chairperson:  J.  Alan  George 

m.  Robert  Koury,  David  F.  McAllister,  and  William  J. 

Stewart*,  "Block  Iterative  Methods  with  Aggregation  for 
Solving  Nearly  Completely  Decomposable  Markov  Chains." 

P.  M.  Dearing,  D.  R.  Shier,  and  D.  D.  Warner*,  "Maximal 
Chordal  Subgraphs  and  Derived  Matrix  Splittings." 

Linda  Kaufman,  “Fast  Direct  Methods  for  Other  Sparse 
Matrix  Problems." 


3B.  Contributed  Papers.  Druid-Dorchester  Room 
Chairperson:  Carlos  A.  Felippa 

m.  Ken  Kanexo,  "The  Turn-Back  LU  Procedure  for  Computing  a 
Sparse  and  Banded  Basis  of  the  Null  Space." 

M.  T.  Heath,  R.  J.  Plemmons*,  and  R.  C.  Ward,  "Sparse 
Orthogonal  Schemes  for  Structural  Optimization  'Jsing  th 
Force  Method." 

A.  S.  Rao,  “Sparse  Matrices  in  Simplification  ot  Systems. 


1:10 


Lunch,  inforiul  afternoon  sessions,  recreation. 


Monday,  October  -’5** 

Session  4A.  Contributed  Papers.  St.  George  Room 
Chairperson:  Conn  : 

7:30  p.m.  Tom  Manteuffel*  and  Vance  r‘aber',  "Conjugate  Gradient 
Imp  1 i es  Norma ! . " 

7:50  Anne  Greenbaam,  Analysis  of  a  Multi-grid  Method  as  an 

Iterative  Technique  for  Soi  ving  Linear  Systems.1, 

3:10  A.  Brandt,  S.  McCormick*,  and  ).  Rune,  "Algebraic 

Multigrid  (AMG)  Tor  Geodetic  conations." 

Session  4B.  Contributed  Papers.  Druid-Dorcnester  Room 
Chairperson:  Jack  J.  Dongarra 

7:30  p.m.  D.  A.  Calahan,  "Sparse  Direct  Methods  for  Vector 
Multiprocessors. " 

7:50  Ronald  D.  Coleman*  and  Edward  J.  Kushner,  "The  Sparse 

Matrix  Library  for  the  FPS-164  Attached  Processor." 

8:10  A.  Sameh*  and  C.  Taft,  "Preconditioning  Strategies  for  the 

Conjugate  Gradient  Algorithm  on  Multiprocessors." 


Break . 


Session  5.  Invited  °aper.  St.  George  Room 
Chairperson:  David  S.  Scott 

3:45  p.m.  B.  N.  Parlett ,  "Software  for  Sparse  Eigenvalue  Problem 


End  of  evening  so.^s’  . 
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Session  6. 


8:45  a.m. 


10:15 


Tuesday,  October  26** 

Invited  Papers.  St.  George  Room 
Chairperson:  Robert  J.  Plemmons 

Michael  T.  Heath,  "Numerical  Methods  for  Large  Sparse 
Linear  Least  Squares  Problems." 

Allen  J.  Pope,  "Geodetic  Computations  and  Sparsity." 


Coffee  break. 


Session  7A.  Contributed  Papers.  St.  George  Room 
Chairperson:  William  G.  Poole,  Jr. 

10:40  a.m.  Aram  K.  Kevorkian*,  Fred  G.  Gustavson,  and  Gary  D. 

Hachtel ,  "A  New  Theory  on  Permuting  Matrices  to 
Block  Triangular  Form." 


11:00 


11:20 


11:40 


Fred  G.  Gustavson*,  Gary  D.  Hachtel,  and  Aram  K . 
Kevorkian,  "An  Improved  Assignment  Algorithm." 

Gary  D.  Hachtel*,  Fred  G.  Gustavson,  and  Aram  K. 

Kevorkian,  "A  1-Pass  Procedure  for  Maximal  Assignment 
and  Finding  the  BLT  Form." 

Fred  G.  Gustavson,  "A  New  Version  of  Tarjan's  Strong 
Connect  Algorithm." 


Session  7B.  Contributed  Papers.  Druid-Dorchester  Room 
Chairperson:  Gene  H.  Golub 

10:40  a.m.  David  S.  Scott,  "LAS02-Sparse  Symmetric  Eigenvalue  Package 


11:00 


11:20 


11:40 


Jane  Cullum*  and  Ralph  A.  Willoughby,  "An  Accelerated 
Block  Lanczos  Procedure  for  Extreme  Eigenvalues  of 
Symmetric  Matrices." 

Paul  o.  Jensen,  "A  Production  Ei genanalysi s  System  for 
Large,  Generalized  Symmetric  Problems." 

Horst  D.  Simon,  "The  Lanczos  Algorithm  with  Partial 
Reorthogonal ization  for  the  Solution  of  Nonsymmetric 
Linear  Systems." 


Tuesday,  decoder  20** 


Session  8A.  Contributed  Papers.  st.  George  Poor) 

Chairperson:  William  G.  Poole,  Jr. 

12:10  p.m.  Silvio  Ursic,  "Inverse  Matrix  Kepresentat i on  with  One 
Triangular  Array  (Implicit  Gauss)." 


12:30 


12:50 


John  de  Pill  is*  and  Wilhelm  Ni ethammer,  "A  Stationary 
Iterative  Method  that  Works  Even  for  Hermitian  Ax  -  b. 

C.  R.  Johnson,  "The  Classes  Mn  ^  and  Properties  of  Spars 

Matrices  Resembling  Those  of  Matrices  in  a  Smaller 
Dimens i on . " 


Session  8B.  Contributed  Papers.  Dru i d-Dorchester  Room 
Chairperson:  Gene  H.  Golub 

12:10  p.m.  Alan  George  and  Joseph  Liu*,  "Row  Ordering  Schemes  for 
Sparse  Givens  Transformations." 


12:30 


12:50 


Alan  George  and  Esmond  Ny*,  "Solution  of  Sparse  Under- 
determined  Systems  of  Linear  Equations." 

M.  S.  Kamel*  and  K.  Sinynal,  "An  Efficient  Algorithm  for 
the  Factor i zat i on  of  Possibly  Rank  Deficient  Matrices. 


Lunch,  informal  afternoon  sessions,  recreation. 


Session  9A.  Contributed  Papers.  St.  George  Poor, 

Chairperson:  Linda  Kaufman 

7:30  p.m.  Iain  S.  Duff  and  John  -;.mc!*,  'The  Multifrontal  Sol.iti 
of  Un symmetric  Sets  of  Linear  Equations. " 

7:50  John  R.  Gilbert*  and  o r; * > t  Schreiner,  "Nested  Dissect'0 

with  Partial  Pivoting. 

8:10  Albert  M.  Erisman.  "Mat.-'  «  Modi  f  i  tat  ion  and  Part  if  ion ’no 


Tuesday,  October  26** 


Session  9B.  Contributed  Papers.  Dru ld-Drochester  Room 
Chairperson:  Yueh-er  Kuo 

7 : 30  p.m.  A.  Behie  and  P.  A.  Forsyth  Jr.*,  "Incomplete  Factori  zat i on 
Methods  tor  Fully  Implicit  Simulation  of  Enhanced  Oil 
Recovery . " 

/ : SO  W.  I.  Bertiger*,  D.  A.  Calahan,  P.  T.  Woo,  "The  Effect  of 

Computer  Architecture  on  Direct  Sparse  Matrix  Routines 
in  Petroleum  Reservoir  Simulation." 

3:10  W.  P.  Kamp,  "A  Three  Dimensional  Finite  Element  Magneto- 

telluric  Modeling  Program." 


8:30  Break. 


Session  10.  Invited  Paper.  St.  George  Room 
Chairperson:  Jane  K.  Cullum 

8:45  p.m.  W.  M.  Coughran,  Jr.,  W.  Fichtner,  E.  H.  Grosse,  and 

D.  J.  Rose*,  "Numerical  Simulation  of  VLSI  Circuits." 


9:30  End  of  Evening  Session. 


Wednesday,  October  27 


Session  11.  Invited  Papers.  St.  George  Room 
Chairperson:  Robert  E.  Funderlic 

3:45  a.m.  Dianne  P.  O'Leary,  "Solving  Mesh  Problems  on  Parallel 
Processors . " 

9:30  Philip  E.  mil,  Walter  Murray,  Michael  A.  Saunders*,  and 

Margaret  ti.  Wrignt,  "Sparse  Matrix  Methods  in 
Opt  ini zat i on . " 


10:15 


Coffee  break. 


Session  12A.  Contributed  Pa^t-r  . 


10:40  a 

11:00 

11:20 

11:40 

Session 

10:40  a 

11:00 

11:20 

11:40 

12:00 

Session 

12:10  p 

12:30 

12:50 


Chai  rperson :  a.  - 

m.  Batman  Nour-Jmij*  i - .  :  ---m"  .  >3  r  -  rn  •  ■  - 

tionin  j  for  Ndl.tmm  oi  •  -  •  *  '  'e.vn*  :  m'  :  s  . 

A.  Jennings*  and  M.  A.  ! •  - 1 e  '•‘etn  ~>d  , 

Solving  AA<  -  b. 

Tony  F.  Chan*  and  kennetn  .a->  w .  Nm;  i  <  n,.dr!y 
Preconditioned  -rvlo,  s  ms  \k>-  Vet»vj  is  tor 
Discrete-Newton  M  .on  •  •• 

Batman  Nour- jmi  d*  and  2erss‘.'.i  ;  jrlii'  ,  'L  1  e-'v.t 
Pre-cond  1 1  i  or  i  ng  . 11 


12B.  Contributed  Papers.  Dm,  •  ouster  Room 

Chairperson:  Margaret  n.  *ri  go* 

m.  R.  J.  Hanson  and  K.  L.  :iiebert*.  A  Snarse  linear 
Programmi ng  Subprogra n. 

Walter  Murray,  "  Null -Spare  Met  nods  ‘or  Large-Scale 
Quadratic  Programing.' 

S.  Thomas  McCormicK ,  "A  Fa«-t  Algorithm  That  Maices  Matrices 
Optimally  Sparse." 

Michael  Engquist,  'A  Par:  :  p  i  o<-  r  ;  Approach  t  o  Pr  x-essiM, 
Networks  and  the  Sol jt ion  of  . parse  Systems. 


Break . 


13A.  Conti'  i  bated  Pa;  -r  - .  *  . 

Chairperson:  "  ber:  M. 

m.  Larry  F.  Bennett,  m  ■  .*.••••  ;*  '  .-•>  ■  e:.t 

Methods  and  >•  ■  •  .>»•  ■■  ‘  ■-  ■'  1  ■  m  va 1 1  : 

Equations  and  l.e  •.  -.t  s  . 

Daniel  B.  Szyld*  -ms  ■  .s*/,  "  i 

on  Sparse  Mu t  r  i  •  ■■  , . 

Jian-Xm  Dec,  >r 

Genera  i  i  zed  -  :  i- ’ 


Wednesday,  October  27 


Session  138.  Contributed  Papers.  Oru i d-Dorchester  Room 
Chairperson:  Margaret  H.  Wright 

12:10  p.m.  Thomas  F.  Coleman,  "Software  for  Sparse  Matrix 
Estimation." 

12:30  Paul  H.  Calamai*  and  Andrew  R.  Conn,  "Continuous  j?.p  Norm 

Location  Problems:  A  Second-Order  Approach." 

12:50  Noel  E.  Cortey,  "A  New  Method  for  Solving  Systems  of 

Linear  Inequalities." 


1:10  End  of  Symposium. 


Poster  Session.  Wi 1 shi re-Canterbury-Wi ndsor  Room 

(Posters  on  display  Monday  and  Tuesday) 


D.  A.  Calahan,  "Performance  of  Sparse  Equation  Solvers  on  the 
CRAY-1." 


Iain  S.  Duff,  Roger  G.  Grimes,  John  G.  Lewis,  and  W.  G.  Poole, 
Jr.,  "Sparse  Matrix  Test  Problems." 


R.  E.  Funderlic  and  R.  J.  Plemmons,  "An  Incomplete  Factori zat i on 
Method  for  Singular  Irreducible  M-Matrices." 


James  E.  Giles,  "Implementation  of  a  Large-Scale  Linear 
Programming  Package  on  a  Minicomputer." 


Roger  G.  Grimes,  John  G.  Lewis,  and  William  G.  Poole,  Jr., 

"Program  for  the  Comparison  of  Reordering  Algorithms  for  the 
Solution  of  Unsymmetric  Sparse  Systems  of  Equations." 


Fred  G.  Gustavson,  "An  Algorithm  for  Computing  a  Full  Assignment 
for  A  Sparse  Matri x. " 


William  A.  l.oden,  'A  Comparison  of  Two  Sparse  Matrix  Processing 
Techniques  for  Structural  Analysis  Applications." 


Monday,  October  25 


St.  George  Room  Druid-Dorchester  Room  Wi 1 shi re-Canterbury-Windsor  Room 


Couyhran,  W.  Fitchner 
Urusse,  D.  J.  Rose* 


Wednesday,  October  27 
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INVITED  ABSTRACTS 


NUMERICAL  SIMULATION  OF  VLSI  CIRCUITS 


William  M.  Coughran,  Jr. 
Wolfgang  Fichtner 
Eric  H.  Grosse 
Donald  J.  Rose* 

Bell  Laboratories 
Murray  Hill,  New  Jersey  07974 


Circuit  analysis  has  motivated  interesting  numerical  analysis  for 
several  decades.  We  will  review  previous  work  with  particular  emphasis 
on  the  type  and  structure  of  the  nonlinear  and  linear  equations  that 
arise.  We  then  suggest  a  new  formulation  of  the  circuit  equations 
based  on  a  hierarchical  structuring  of  the  circuit.  Our  reformulation 
of  the  equations  leads  naturally  to  the  circuit  substructuri ng  that  is 
often  termed  "macromodeling",  currently  in  vogue. 

We  will  present  several  numerical  examples  of  tnacromodel i ng.  In 
particular  we  will  emphasize  the  construction  of  a  transistor  device 
macromodel  directly  from  a  device  simulator  as  well  as  from  the  more 
traditional  equivalent  circuit. 


Session  10.  8:45  -  9.30  p.m, ,  Tuesday 


♦Indicates  the  presenter  of  the  paper 


than  one  aathor  is  given, 


DIRECT  METHODS  FOR  SOLVING 
SPARSE  SYSTEMS  OF  LINEAR  EQUATIONS 


Iain  S.  Duff* 

John  K.  Reid 
AERE  Harwel 1 
Oxon,  0X11  ORA,  ENGLAND 


We  survey  algorithms  and  software  for  solving  sparse  systems  of 
linear  equations  paying  particular  attention  to  recent  developments. 
We  classify  the  various  algorithms  according  to  the  type  of  system  they 
solve  (i.e.,  unsymmetric,  symmetric  definite,  symmetric  indefinite, 
unsymmetnc  but  with  symmetric  pattern)  and  whether  they  perform 
pivoting  for  numerical  stability.  We  consider  >oth  algorithms  which 
work  in  main  memory  and  those  which  use  auxiliary  storage. 

We  illustrate  the  performance  of  most  of  the  software  we  discuss 
by  runs  on  test  problems  and  give  a  critical  evaluation  of  each, 
stressing  their  strengths,  weaknesses  and  restrictions. 


Session  1.  9:00  -  9:45  a.m.,  Monday 


ITERATIVE  METHODS  FOR  SOLVING 
LARGE  SPARSE  LINEAR  SYSTEMS 


Std.fi  t-y  C.  .  •’  sen-'.t 
Y  d  1  r1  '  n  l  v  5  r  S  i  t.  y 

New  Haven.  Connecticut  on  5  20 


This  talk  will  (attempt  to)  surve.  the  state  jf  the  art 

iterative  and  semi -i terat  i  ve  methods  *or  solving  large  sparse  systems 

of  linear  equations  Ax  =  b.  Tnree  classes  of  problems  will  be 
considered,  corresponding  to  the  coefficient  natn,  being  symmetric 
positive  definite,  symmetric  indefinite,  or  nonsymmet ri c.  For  each 

class,  the  emphasis  will  be  on: 

•  Krylov  subspace  methods  -  which  choose  the  mth  iterate  from  the 

Krylov  subspace  Span  lb,  Ab,  A?o,  Ain_  1  b  ’•  and  thus  only  require 

the  ability  to  multiply  a  vector  by  the  coefficient  matrix  (e.g., 

the  conjugate  gradient  and  Chebyshev  semi - i terat i ve  methods); 

•  preconditioning  techniques  -  netnods  for  ’rescaling"  the  original 

system  to  increase  the  rate  of  convergence  (e.g.,  incomplete 
factorization); 

•  reduction  techniques  -  methods  for  reducing  a  problem  m  one 

class  to  a  smaller  problem  in  the  same  class  (e.g.,  cyclic 

reduction)  or  a  problem  in  another  class  (e.g.,  forming  the  normal 
equations)  while  maintaining  whatever  sparseness,  symmetry,  and 
definiteness  was  present  in  the  original  system; 

•  general-purpose  software. 

The  written  version  of  the  talk  will  include  an  extensive  bibliography 
of  recent  papers. 


Session  1. 


9:45  -  10:30  a.m.,  Monday 


SPARSE  MATRIX  METHODS  IN  OPTIMIZATION 


Phi  lip  E.  Si  1 1 
Walter  Murray 
Michael  A.  Saunders* 
Margaret  H.  Wright 
Stanford  University 
Stanford,  California  94305 


Optimization  algorithms  typically  require  the  solution  of  many 
systems  of  linear  equations  Bkyp  -  bk.  When  large  numbers  of  variables 
or  constraints  are  present,  these  linear  systems  could  account  for  much 
of  the  total  computation  time. 

Both  direct  and  iterative  equation  solvers  are  needed  in 
practice.  Unfortunately,  most  off-the-shelf  solvers  are  designed  for 
single  systems,  whereas  optimization  problems  give  rise  to  hundreds  or 
thousands  of  systems.  To  avoid  ref  actor i zi ng ,  or  to  speed  the 

convergence  of  an  iterative  method,  it  is  essential  to  note  that  Bk  is 
related  to  Bk_i . 

We  review  various  sparse  matrices  that  arise  in  optimization,  and 
discuss  compromises  that  are  currently  being  made  in  dealing  with 
them.  Since  significant  advances  continue  to  be  mad'  ith 
system  solvers,  we  give  special  attention  to  methods  :t,..  allow  b^ch 
solvers  to  be  used  directly  on  modified  systems  (e.g.,  the  PF]  update; 
use  of  the  Schur  complement).  At  the  same  time,  we  hope  that  future 
improvements  to  I i near-equati on  software  will  be  oriented  more 
specifically  to  the  case  of  related  matrices  . 


Session  11.  9:30  -  10:15  a.m.,  Wednesday 


NUMERICAL  METHODS  TOR  LARGE  SPARSt 
LINEAR  LEAST  SQUARES  PROBLEMS 


M  i  c  n  a  e !  T .  1  •  a  •_  n 

Union  Carbide  Corporation,  Vjc  ’  ear  Pi;;-.’ 
Oak  P i doe ,  Tennessee 


La  rye  sparse  least  squares  problems  arise  in  .nany  appl i cat i r  s , 
including  geodetic  network  adjustments  and  finite  element  strict  mjl 
analysis.  Although  geodesists  and  engineers  have  been  solving  Sjci 
problems  for  years,  it  is  only  relatively  recently  that  numerical 
analysts  have  turned  attention  to  them.  In  this  talk  we  present  a 
survey  of  numerical  methods  for  large  sparse  linear  least  squares 
problems,  focusing  mainly  on  developments  since  the  last  comprehensive 
surveys  of  the  subject  published  in  1976.  We  consider  direct  methods 
based  on  elimination  and  on  orthogona  1  i  zat  i  on ,  as  well  as  various 
iterative  methods.  The  rami f i cat i ons  of  rank  deficiency,  constraints, 
and  updating  are  also  discussed. 

Session  6.  8:45  -  9:30  a.m.,  Tuesday 


SOLVING  MESH  PROBLEMS  ON  PARALLH  PROCESSORS 


'Marne  ejr, 

!■  :  vers  *  ty  o*  Mar  /  ’  a  no 
;e  Park,  M<jr/ :  ant! 


•  i;  :-»s*»  t n  a  t  tr-e  sparsity  structure  ot  an  opti.ii  zation  problem  or 
ot  linear  ,jr  non’ inear  equations  corresponds  to  a  rectangular 
,ru  ot  nodes  ’n  whi^n  each  node  is  linked  to  any  or  all  of  its  north, 
so'j'.n,  east,  west,  northeast,  northwest,  southeast,  and  southwest 
n e i  i n o o r s ,  Seen  problems  arise  in  fields  including  discretization  of 

aartial  ni f f erent i a i  equations,  network  problems,  and  image 

processing.  I*  the  problem  is  to  be  solved  on  a  parallel  computer 
usmg  an  iterative  method,  the  Jacobi  method  allows  parallelism,  but 
special  orderings  of  nodes  are  required  to  exploit  parallelism  in 
methods  such  as  Gauss-Sei del ,  SOR ,  or  conjugate  gradients. 

In  this  worx,  we  discuss  requirements  on  parallel  computer 
architectures  and  algorithms  which  permit  efficient  solution  of  large 
sparse  problems.  Node  orderings  and  processor  arrangements  are 
presented  which  assign  each  processor  to  a  small  number  of  points  and 
enable  each  to  work  in  parallel  with  other  processors  and  only  limited 
data  transfer  among  processors.  Examples  discussed  include  the  nine- 
point  finite  difference  operator  and  an  optimization  problem  with  bound 
constraints. 


Session  11.  8:45  -  9:30  a.m.,  Wednesday 


SOFTWARE  FOR  SPARSE  EIGENVALUE  PROBLEMS 


Berestord  N.  Parlett 
University  of  California 
Berkeley,  California  14  : 


There  are  several  good  packages  for  solving  sparse  linear  systems, 
yet  nothing  comparable  for  eigenvalue  calculations.  Why  not? 

This  question  will  be  discussed  along  with  a  survey  of  what  is 
available  to  the  public.  Some  good  software  is  buried  in  applications 
packages  in  engineering,  chemistry,  and  physics  centers. 

Some  judgments  and  comparisons  will  be  offered  concerning  subspace 
iteration  and  versions  of  the  Lanczos  algorithm.  Of  more  importance  is 
improved  understanding,  and  consensus,  among  the  experts  since  1973. 

It  seems  as  though  eigenvalue  programs  will  have  to  be  specially 
crafted  for  computers  such  as  the  Cray  1.  This  prospect  poses  the 
question  of  the  right  level  of  portability  to  contemplate  for  software 
which  is  going  to  make  more  demands  on  the  host  system  than  does  any 
subroutine  in  EISPACK.  The  use  of  secondary  storage  is  a  key  factor  in 
the  efficiency  of  a  program  and  this  makes  it  difficult  to  write 
software  which  is  both  effective  and  independent  of  the  operating 
system. 


Session  5.  8:45  -  9:30  p.m. ,  Monday 


I 


GEODETIC  COMPUTATIONS  AND  SPARSITY 


i 


Alien  J.  Pone 
National  Geodetic  Survey 
NOAA/NOS 

Rockville,  Maryland  20352 


The  least  squares  adjustment  of  geodetic  networks  is  a  natural 
arena  for  the  application  of  sparse  matrix  technology.  In  this  taU, 
geodetic  computations  are  reviewed  for  the  non-geodesist.  Special 
emphasis  is  given  to  those  features  (of  which  there  are  many)  that  set 
these  computations  apart  from  most  other  sparse  problems.  A  survey  of 
the  geodetic  involvement  with  sparsity  begins  near  the  origins  of  the 
science  and  continues  through  recent  trials  of  various  re-ordering 
strategies.  Although  it  can  be  maintained  that  the  exploitation  of 
sparsity  has  reached  a  certain  palteau,  there  still  remain  numerous 
challenges.  In  particular  these  arise  from  the  large  size  of  some 
geodetic  data  sets,  such  as  that  involved  in  the  re-adjustment  of  the 
North  American  datum. 


Session  6.  9:30  -  10:15  a.m.,  Tuesday 


►V 


CONTRIBUTED  ABSTRACTS 


INCOMPLETE  FACTORIZATION  METHODS  IOK  TUNY 
IMPLICIT  SIMULATION  OF  ENHANCED  OIL  RFCOVi  RY 


A.  ben  10 

P.  A.  Forsyth  Jr.* 
Computer  Modelling  iroup 
Calgary,  Alberta,  Canada  T/l.  /Ad 


Fully  implicit  simulation  of  enhanced  oil  recovery  using  thcma! 
(steam  and  in  situ  combustion)  .nethods  gives  rise  to  a  highly  struc¬ 
tured  block-banded  non -symmetric  Jacobian.  Solution  of  large  pron’ens 
requires  effective  iterative  methods. 

Recently,  incomplete  factori  zat  1  on  nethods  (ILU)  have  been  use  i  t.i 
solve  these  systems.  In  this  paper,  several  ILU  methods  are  developed 
using  natural,  diagonal  (D2)  and  alternating  diagonal  (04)  ordering. 
Diagonal  ordering  was  first  used  for  symmetric  systems  by  Watts,  while 
Tan  and  Letkeman  suggested  alternate  diagonal  ordering  for  non- 
symmetric  problems.  The  method  of  Watts  is  generalized  to  the  strongly 
non-symmetric  case,  and  an  improvement  to  the  algorithm  of  Tan  and 
Letkeman  is  developed  which  saves  a  considerable  amount  of  wo rk.  In 
addition,  several  degrees  of  f actor i zat i on  are  used  for  all  these 
orderings.  These  techniques  are  all  accelerated  with  the  URTH.OMIN 
al gori tnm. 

F.ach  of  these  methods  is  developed  with  vector  machines  in  mind, 
and  particular  attention  is  paid  to  those  portions  of  the  algorithm 
which  can  be  readi  ly  vectorized.  All  the  ILU  methods  can  also  be  used 
with  the  COMBINATIVE  technique. 

Results  are  presented  for  several  model  problems  and  test 
Jacobians  generated  from  steam  simulations.  The  results  snow  tne 
effects  of: 
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•  ■  ,nt  er  ■>  p -  ‘  •  (t>  Iterative  ‘•S't.  n  )  1  >)  f  >r  r.ht.*  sol  it  i  on  o*  n  -  n  natrix 

-Mat;')'",  ar  !  i  n>\t  r  least  s  r.jares  profile  ss  ,  and  accelerated  general¬ 
ize.]  row  project  ion  nefuids  pert  enrect’ve  Iterative  Methods)  for  the 
so!  :.i  1 1  o r  of  ’!  -  n  "ia  t  r  i «.  .•>  p;a*  ions  and  linear  least  squares  problems. 
-! !  project’,  ve  net  nods  discussed  are  finite  -net  nods,  so  that  line  the 
net  nod  o“  cone. gate  gradients,  they  are  guaranteed  to  converge  to  the 
exact  sol ut ion  sought  ir  a  ‘inite  number  of  steps  i*  infinite  precision 
arithmetic  were  possible.  Algorithms  and  software  for  computerized 
implementation  of  the  methods  in  solving  sparse  matrix  equations  and 
least  squares  problems  as  well  as  operation  counts  are  mentioned. 

Important  properties  of  the  new  methods,  such  as  those  which  follow  are 
stressed.  The  methods  appear  to  work  comparatively  well  on  any  matrix 
least  squares  problem  or  matrix  equation  of  the  form  Ax  =  h,  where  A  is 
any  complex  or  real-valued  m  <  n  matrix  and  b  is  a  complex  or  real¬ 
valued  m-tuple.  No  special  hypothesis  on  the  structure  of  the  matrix  A 

is  needed.  For  large  seal-*  natrix  problems,  convergence  to  an  accept¬ 
able  solution  often  occurs  ouch  sooner  than  predicted  for  the  exact 
solution.  The  methods  can  produce  extremely  accurate  results,  and  in 

the  case  of  least  squares  problems,  do  not  appear  to  increase  ill -con¬ 
ditioning  as  normal  equations  sometimes  do.  The  methods  may  be  further 
accelerated  in  some  cases  by  using  relaxation  factors.  In  many  cases, 
single  precision  u^ithuet'c  may  be  employed  to  generate  accurate 
solutions,  even  in  the  case  of  ill-conditioned  systems.  The  rate  of 
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THE  EFFECT  OF  COMPUTER  ARCHITECTURE  ON  DIRECT  SPARSE 
MATRIX  ROUTINES  IN  PETROLEUM  RESERVOIR  SIMULATION 


W.  I.  Hertiqer* 

Chevron  Oi  Field  Researc'i  Company 
La  Hahra,  California  A06ol 

D.  A.  Cal  aha n 
University  of  Mi  chi q an 
Ann  Arbor,  Michigan  48109 

P.  T.  Woo 

Chevron  Oil  Fie'd  Pose-arch  Company 
La  Habra,  California  9df>31 


Large  systems  of  linear  equations  are  solved  daily  in  reservoir 
simulation.  These  equations  arise  from  finite  difference  approxima¬ 
tions  to  systems  of  partial  differential  equations  describing  the  flow 
of  fluids  in  a  reservoir.  With  the  advent  of  new  vector  computer  tech¬ 
nology  the  algorithms  for  solving  the  linear  equations  must  be  changed 
to  make  better  use  of  the  machine  archi tecture. 

In  solving  linear  equations  on  a  vector  computer,  there  is  a  trade 
off  between  computation  rate  and  work  (the  amount  of  arithmetic).  At 
this  stage  of  development,  optimum  use  of  the  vector  computer  has  led 
us  to  change  from  a  Vale  general  purpose  sparse  matrix  routine  with  the 
inner  loops  coded  in  assembly  language  to  an  assembly  coded  bandsolver 
developed  at  the  University  of  Michigan.  The  sparse  routine  is 
optimized  to  do  the  least  amount  of  arithmetic  and  is  best  suited  for  a 
scalar  computer.  The  bandsolver  is,  on  the  other  hand,  optimized  to 
achieve  maximum  computation  rate  on  a  vector  computer.  The  trade  off 
for  solving  linear  equations  in  reservoir  simulation  is  in  favor  of  the 
bandsolver. 

When  equations  involving  the  wellbore  pressures  are  added  to  the 
flow  equations,  the  added  portion  of  the  matrix  is  not  banded.  We  will 
describe  how  the  bandsolver  can  still  be  us*-!  effectively  with  an  algo¬ 
rithm  attributed  to  George.  Ci  rr'ivirs  ‘--e  tween  the  sparse  .nd  pjr  j- 
solvers  based  on  operation  counts  and  --ct.ia  i  timing  rmb  nr  the  Cray  - 1  :> 
will  be  presented  for  model  and  real  reservoir  problems.  ;<esu 1 1$  f  or 
the  Cray  XMP  wi  1  I  be  presented  ’f  the  oeoenmar*  r/o  can  oe  /-..imp  1  eto-i 
in  time.  Future  improvements  to  roe  direct  sol  ,f'or>  ‘  *' 

equations  on  vector  computers  will  :>e  Hi  •>  cn  ■■■  w  -i . 
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ALGEBRAIC  MULJIGRIl)  (AM6)  FOR  GIODLTIC  EQUATIONS 
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i-ort  ;o  i  1 1  ms ,  Colorado 


Mu  1  ti gri a  a  1  joritnm  do  si  gn  regu  i r  es  that  its  major  processes  be 
tailored  to  each  app  1 1  cat  i  o” .  "his  can  otter  be  done  manual 1y  at  the 
tHe  the  di  scret  i  zd  t  i  on  cocertore  is  r>t»ing  designed,  but  not  always. 
The  engineer  'nay,  for  example ,  preselect  a  fine  arid  finite  element 
tri  angulation  so  that  detenu  m  oy  even  what  the  coarser  grids  should 
be  is  difficult  if  not  impossible  .  Moreover,  some  problems,  such  as 

those  that  arise  in  geodetic  surveying,  are  inherently  discrete.  So 
it  is  important  to  consider  methods  for  automatic  mul  tigrid  design. 

Algebraic  mul  1 1  grid  (AM'S)  is  one  such  method,  where  the  discrete 
problem,  represented  either  a  matrix  or  operator  stencil,  is 

subjected  to  a  preprocessing  stage.  A MG  bases  its  decisions  on  the 
concept  of  strong  coupling  whicn,  loosely  speaf-ing,  is  a  way  to 
interpret  effective  dependencies  between  variables  via  coefficients 

in  the  discrete  operator  tnat  connect  them.  Tens  concept  provides  a 
motivation  for  deciding  whicn  variables  are  to  constitute  the  coarse 
“grid"  and  for  determining  interpolation.  The  remaining  multiynd 
processes  can  be  developed  from  this  design  in  a  fairly  straiyht- 

f  o  rw  a  r d  way  . 

Tms  is  a  pre!  i:r.j  nary  report  on  experience  with  several  versions 

of  A  MG  applied  to  el  1 ipti-  differential  problems  that  include  diffusion 
equations  (with  coefficients  that  are  both  widely  varying  and  strongly 
discontinuous) ,  anisotroo’  fluid  "lows.  and  non-uniform  grids. 
Special  attention  :  s  ;i  to  the  ;.ntentia!  ot  •'•'•I.,  tor  solving  purely 
algebraic  , mo, v-  "r,.  ]  ••  ; .  •  geod-M  c  editions. 
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SPARSE  DIRECT  METHODS  FOR  VECTOR  MULTIPROCESSORS 


0.  A.  Calahan 
University  of  Michigan 
Ann  Arbor,  Michigan  4,-5109 


The  CRAY  X-MP  is  a  harbinger  of  a  class  of  forthcoming  AljAFLOP 
supercomputers  composed  of  arrays  of  vector  processors  connected  to  a 
common  main  memory.  Related  sparse  matrix  Methods  must  combine  the 
parallel  (multiprocessor )  and  vector  taxonomies  that  nave  been  studied 
independently  in  the  past. 

Computationally,  the  solution  may  be  phrased  as  a  two-level 
synchronization. 

1.  The  vectorized  inner  loop,  precisely  synchronized  at  the 
instruction  level,  and  operating  independently  among  the 
processors . 

2.  The  parallel  outer  loop,  approximately  synchronized  at  the 
instruction  block  level. 

If  the  vectors  are  sufficiently  long,  both  loops  may  be  precisely 

synchronized  by  distributing  the  inner  loop  vectors  among  the 
processors  (outer  loop)  and  then  having  them  operate  in  an  S I MD  mode. 
This  is  undoubtedly  the  manner  in  which  large  full  and  banded  systems 
will  be  solved,  but  is  al yorithmical ly  uninteresting. 

Rather,  consider  the  case  of  a  randomly-sparse  matrix  where  each 
element  is  a  rectangular  block  representing  coupling  of  unknowns  at  two 
grid  points  in  a  finite  element  or  finite  difference  representation. 
Further,  these  blocks  are  different  sizes  corresponding  to  the 

different  number  of  unknowns  at  each  grid  point,  and  may  change 

dynamically  in  a  time-dependent  problem. 

It  is  proposed  that  a  solution  proceed  by 

1.  Local  decoupling  of  olocxs  by  reordering  of  elements 

(consistent,  for  example,  with  nested  dissection ; ,  me 

2.  Dynamic  scheduling  of  independent,  block  processing  m  ng  toe 
processors;  dynajnic  scned.i  t  1  rg  -s  necessary  to  com:  o-.af-'  r  .r 
different  blocx  sizes. 

An  extended  CRAY  X-MR  arch  1 1  oct  ,  e-a-med  t  e  im.ot’-:!  ;«►.  > »  us¬ 
ance  of  such  an  algorithm.  A  simulator  ~s  being  developed  'nr  tm ; s 
machine,  and  wi  1 1  quite  likely  permit  some  precise  ore ]ect i ons  by 
time  of  the  conference.  Actual  runs  o'"  a  n-'ssir.le. 
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PERFORMANCE  OK  SPARSE  EQUATION  SOLVERS  ON  THE  CRAY-1 
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A 1  y  o n  i  t  n  t) arid  pert  jr.'yn  .e  o'  live  ;  I  asses  of  sparse  •  ta  t  r  i  x 
solvers  are  di scussed . 

1.  General  sparse  solvers.  These  solve  randonl y-sparse  equations 
described  symbolically  in  ~o  1  umn-ordereJ  for~i  and  are  compatible  with 
traditional  scalar  paoages.  The  approaches  investigated  both  involve 
two-step  symbol  ic/nii-neric  processing  (in  the  manner  ot  Gustavson,  et  al). 

(a)  Hybrid  vector-seal ar  solution  of  medium-density  systems.  The 
symbolic  phase  identities  dense  segments  of  columns  with  more 
than  a  user-prescribed  value  (n)  o*  contiguous  non-zeros. 
These  segments  are  processed  in  vector  mode  and  the  remainder 
in  scalar  mode.  Maximum  rate:  o’  MKlQhb. 

(b)  Decoupled  scalar  sol  .it ion  of  hignly- sparse  systems.  A  varia¬ 
tion  of  code  generation  methods  for  a  pipelined  processor 
involves  tne  local  decoupling  of  hi jh iy-sparse  equations  (in 
the  manner  of  nested  dissection).  Scalar  machine  code  is  then 
generated  which  see* s  to  "cram"  the  floating  point  pipelines 
with  independent  scalar  code.  Tyi)ical  rate:  lb  MFLOPS. 

2.  Special  sparse  solvers.  Higher  performance  can  be  achieved 
when  special  sparsity  structure  can  be  identified  by  the  user.  In 
general,  vectors  car  he  -i  s-- f  .<  1 T  v  def  i  ni.a  ei  ?  »->er  within  a  dorse  submatr’  x 
or  across  si:ni  lari  /-struct  ,r(-;i  -yih-tutri  cos. 
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CONTINUOUS  J(p  NORM  LOCAT  ION  PROBLLMS:  A  SF  COND-ORDER  APPROACH 


d  ■.  I  H .  ;J  :• 
Andre*  v . 
university  of 

Waterloo,  '■*  or  o. 


A  second-order  a  I  gori  thm  for  sol  vi  ng  a  prototype  continuous  mini- 
sum  multi faci  1  ity  location  problem  involving  distances  is  presented. 
This  problem  corresponis  to  a  special  case  of  'the  p-median  problem  ann 
a  continuous  version  of  the  quadratic  assignment  problem. 

We  demonstrate  how  projectors  can  be  used  to  circumvent  tne 
difficulties  associated  with  the  nondi fferentiaoi 1 i ty  of  these 
problems.  Tne  technique,  which  is  an  extension  of  an  earlier 
first-order  method  used  by  the  autnor,  provides  a  unified,  numerically 
sound  and  stable  approach.  The  implementation  is  the  first  to 
efficiently  exploit,  the  special  structure  of  the  problems  under 
consideration.  In  particular,  I)  we  are  able  to  solve  the  linear 
systems  that  arise  in  the  development  by  stable  means  that  do  not 
suffer  from  fill-in,  2)  a  special  linesearch  which  is  based,  in  part, 
on  the  results  presented  in  a  paper  by  Overton  recognizes  the 
possibility  of  first  derivative  discontinuities  and  second  derivative 
unboundedness  along  descent  directions,  and  3)  the  degeneracies  which 
are  an  inherent  character! sti c  of  these  problems  are  dealt  with  using 
simple  perturbati ons. 

Although  details  are  initially  given  for  only  unconstrained  fixed 
ig  norm  problems  we  show  how  the  framework  can  be  readily  extended  to 
mixed  norm  problems  as  well  as  to  constrained  problems. 
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NONLINEARLY  PRt CONDI  HOMED  KRYLOV  SUM SR ALE  HE  I  HODS 
FOR  DISCRETE  -NEWTON  ALGOR  HUMS 


One  promising  idea  that  sore  PecuV  im j-  •  •  ’  ■  ■ :  i»  the  ,ip;>  1  i  - 
cation  of  Krylov  subspace  net.  rods  *  or  so' .mu  ;  :■->  '-.-v-r  .  v  s  t s  Aw  b 
that  arise  in  the  inner  loop  of  a  Newt  on  - 1  ■  •  ?  'tc-rat'-  *  <j  r 

solving  nonlinear  systems  E'x)  ---  >  (with  A  :v  1  ♦.•■>**  . •».  :sd"  or  >'  •  is 
the  use  of  the  directional  dif f erencinu  c  (  v  u ;  -  -  ;  v ;  for  a.»^ro*i  - 
mating  the  matrix-vector  product  Au  in  the  >;ryi  s  ihs.  methods. 
This  requires  only  function  evaluations  and  j  /  v  n  licit  /*.  i  u  a*- ;  on 
and  storage  of  the  Jacobian.  We  are  interested  m  accelerating  toe 
convergence  rate  of  the  inner  loop-  by  preconditioning  while  still 
employing  directional  differencing.  however,  sioue  most  precondition¬ 
ing  techniques  are  derived  from  the  matrix  elements  of  A  expl icitiy,  it 
is  not  obvious  how  to  apply  them  wnen  A  is  not  <■•«.;>!  i  c.  1 1 1  y  available. 
We  have  deri  ved  an  algorithm  for  precond ’  f  i  cn ■  ng  the  irylov  subspace 
method  with  directional  differencing  that  reduces  to  the  SSOR 
preconditioning  of  A  in  the  linear  case.  It.  on\y  r-eq  li  res  evaluating 
the  diagonal  elements  of  tne  .Jacobian,  which  ca"  be  approximated  by 
function  evaluations.  The  overal 1  al gorithm  is  thus  well-suited  for 
large  and  sparse  problems,  especially  when  f. ns. Moo  evaluations  are  not 
too  expensive.  SomeM  cu  1  experiment.  ,  show  :  n  ■>  t  ;  r  *•  so  H  near  precon- 
iitioning  is  as  effective  a  .  in  to-.  !  inear  mr.fti 
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THE  SPARSE  MATRIX  LIBRARY  FOR  THE  FPS-164  A  IT  ACHED  PROCESSOR 


Ronald  Li.  Coi  enw'* 
Edward  J.  kushner 
Floating  Point  Systems ,  !n. 
Portland,  Oregon  y.’LVd 


This  paper  describes  the  sparse  matrix  routines  tor  t'u-  •  •  ■-t \ n : 
Point  Systems  FPS-164  Attached  Processor  and  presents  associated  'i.-nch- 
mark  results.  The  highly  parallel  and  pipelined  arcni  tec  tun.-  o  t  trie 
FPS-164  provides  high  performance  for  vector  and  matrix  computations  at 
relatively  low  costs.  Included  as  part  of  the  FPS-164  Program 
Development  Software  is  the  APMATH64  Math  Library  which  contains  nearly 
500  subroutines  that  are  organized  into  14  sub-libraries.  One 
sub-library  is  the  Sparse  Matrix  Library  which  contains  13  routines 
for  the  solution  of  sparse  linear  systems  by  both  direct  and  iterative 
methods.  The  Sparse  Matrix  Library  also  includes  12  routines  for 

performing  matrix  arithmetic  operations  with  sparse  matrices.  In 
addition  to  these  routines,  the  Advanced  Math  Library  (also  part  of 

APMATH64)  includes  a  sparse  linear  programming  routine  and  a  profile 

oriented  linear  system  solver  routine.  The  library  routines  are 

written  in  APAL64,  the  FPS-164  assembly  language.  Typically,  these 
routines  execute  15  to  25  times  faster  than  the  equivalent  FORTRAN 
routines  on  the  VAX-11/780. 
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SOFTWARE  FOR  SPARSE  MATRIX  ESTIMATION 


Thomas  F.  Coleman 
Cornell  Jn 1  vers i t y 
Ithaca,  New  Yors.  14853 


The  solution  of  sparse  systems  of  nonlinear  equations  or  qiffer- 
ential  equations  usually  requires  the  determinant  ion  of  the  Jacobian 
matrix  of  a  nonlinear  mapping.  Often  it  is  advantageous  to  estimate 
the  Jacobian  matrix  by  function  (or  perhaps  gradient)  differences. 
When  the  Jacobian  (or  perhaps  Hessian)  matrix  has  a  known  sparsity 
structure,  the  number  of  function  evaluations  needed  for  the  estima¬ 
tion  can  be  quite  small  compared  to  the  dimension  of  the  problem. 
Recently,  Coleman  and  More  have  exploited  a  graph  theoretic  view  to 
develop  efficient  algorithms  to  estimate  sparse  Jacobian  and  Hessian 
matrices.  In  this  talk  we  describe  the  resulting  software  that  has 
been  developed,  at  Argonne  National  Laboratory,  for  the  sparse  matrix 
estimation  problem. 
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A  NEW  METHOD  FOR  SOLVING  SYSTEMS  OF  LINEAR  INEQUALITIES 


Noel  E.  Cjrtyy 
Lockheea-Ca'  >fonna  Co  many 
Burbanx,  California  i L i> 2 ET 


A  new  iterative  descent  method  to  perform  least-squares  is  ,r'v- 
sented.  It  is  shown  how  a  local  weighted  least-squares  solution  can  be 
achieved  by  adding  extra  variables  and  extra  equations  to  a  system  of 
polynomial  equations  and  performing  sequential  linear  weighted  least- 
squares.  When  this  new  method  is  applied  to  solving  systems  of  linear 
inequalities  it  always  leads  to  a  solution  when  tnere  is  one.  When 
there  is  no  solution  it  converges  to  an  approximate  least  squares  solu¬ 
tion.  Since  the  convergence  rate  depends  on  the  initial  guess  it  is 
difficult  to  compare  this  method  with  others.  Nevertheless  it  is 
established  that  globally  the  convergence  rate  is  not  an  increasing 
function  of  m  min  (e,v)  where  e  is  the  number  of  equations  and 

inequalities  and  v  is  the  number  of  variables.  Of  course  this  method 
is  a  good  candidate  to  take  advantage  of  existing  sparse  matrix 
algorithms  for  achieving  linear  1  east-squares  solutions. 
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LS0D28:  A  VARIANT  OF  LSODE  FOR  PROBLEMS  HAVING  A 
GENERAL  LARGE  SPARSE  JACOBIAN 


R.  L.  Cox 

Union  Carbide  Corporation,  Nuclear  Division 
Oak  Ridge,  Tennessee  37830 


LS0D28  represents  a  merging  of  LSODE  (a  state-of-the-art  computer 
program  written  by  A.  C.  Hindmarsh  for  numerically  integrating  stiff 
systems  of  first  order  ordinary  differential  equations)  and  MA28  (a 
package  written  by  I.  S.  Duff  for  solving  sets  of  sparse  unsymmetric 
linear  equations)  which  is  designed  to  handle  stiff  problems  in  which 
the  Jacobian  is  a  general  large  sparse  unsymmetric  matrix.  Because 
MA28  performs  LU  decompositions  incorporating  partial  pivoting  for 
numerical  stability,  LS0D28  avoids  the  requirement  that  the  Jacobian  be 
positive  definite  or  diagonally  dominant  as  demanded  by  some  previously 
existing  codes  which  do  not  incorporate  pivoting. 

The  use  of  sparse  matrix  techniques  as  opposed  to  manipulation  of 
the  full  matrix  results  in  great  savings  in  both  computer  storage  and 
execution  time  for  large  problems.  At  the  same  time,  questions  must  be 
faced  which  are  of  little  or  no  concern  in  the  use  of  full  matrix 
algorithms. 

These  issues  include  an  efficient  user-ori ented  method  for 
representing  the  structure  (location  of  nonzeros)  of  the  Jacobian 
matrix,  the  frequency  at  which  the  pivot  strategy  is  to  be  revised,  the 
possibility  of  changes  in  the  structure  of  the  Jacobian  (zero 
derivatives  becoming  nonzero)  during  the  integration,  and  the  efficient 
generation  of  the  Jacobian  when  its  elements  must  be  calculated 
numeri cal ly . 
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AN  ACCELERATED  BLOCK  i.ANCZOS  PROCEDURE  FOR 
EXTREME  EIGENVAULES  OF  SYMMETRIC  MATRICES 


Jane  Cu  i  !■:",* 

RU  1  ph  e .  will  Onqn  hy 
IBM  T.  J.  watson  Research  Center 
York  town  Heights.  New  Yyr»,  1?)S9>5 


Typically  in  an  iterative  block  Lanczos  procedure,  the  user 
specifies  a  priori  tne  amount  of  computer  storage  available  for  tne 
block  computation.  Therefore,  there  is  a  tradeoff  between  the  number 
of  vectors  included  in  the  first  dIock  and  the  number  of  blocks  that 
can  be  generated  on  a  given  iteration.  The  size  of  the  first  block 
determines  the  gap  between  the  desired  eigenvalues  and  the  eigenvalues 
not  being  approximated  by  the  procedure.  The  more  vectors  used  in  the 
first  block  the  larger  these  gaps.  The  number  of  blocks  used  within  a 
given  iteration  determines  the  effective  spread.  Tne  more  blocks  that 
are  used,  the  smaller  the  effective  spread.  The  convergence  rate 
depends  upon  the  ratio  of  these  gaps  and  the  effective  spread. 

Since  the  user  typically  doesn't  know  the  gaps  a  priori,  generally 
speaking  the  best  overall  approach  is  to  try  to  maximize  the  number  of 
blocks  used  on  each  iteration.  The  number  of  blocks  used  in  each 

iteration  in  the  block  procedure  in  Cullum  and  Donath  [1974]  could  vary 
as  the  computations  proceeded.  The  second  block  was  computed  using  a 

modified  Gram-Schmidt  orthogonal i zat i on  with  pivoting  for  size  and  any 
vector  whose  norm  was  less  than  a  given  program-generated  tolerance  was 
dropped  from  the  block.  This  reduction  in  size  of  the  second  block 

could  lead  to  an  increase  in  the  number  of  blocks  generated  on 

subsequent  iterations  and  an  acceleration  of  convergence.  Such  an 
acceleration  of  convergence  however,  might  not  be  observed  until  after 

a  large  number  of  iterations  and  in  fact  might  never  happen  due  to  the 
particular  sizes  of  the  b 1 oc k s  involved. 

In  this  talk  we  present  an  iterative  block  Lanczos  procedure  that 
maximizes  the  number  of  blocks  generated  on  each  iteration.  The  full 

unnormalized  second  block  is  computed  on  t-acn  iteration.  Trier,  however, 
using  a  selection  procedure  based  upc"  an  opt.  i  r.i  Nation  argument,  or,  1  v 
the  'best'  vector  in  the  second  bloc*  is  a !  1  >*i-i  to  perpet  sate  Use’.*-. 
Succeeding  blocks  are  reortnogona  I  f;*e*i  w.r.t.  tne  ancest  ars  ot  tne 
vectors  dropped  from  trie  second  bi  •>..». .  in  many  this  j.gir,;a-.h 

yields  greatly  accelerated  converge"  o  ::r  what  obtained  b, 

approach  in  Cullum  and  Donat "  ,  1  l  j  .  Curs'  •  c  •  ex.  i  ■  g  ■  i  ■  ■ 

incl uded. 
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MAXIMAL  CHORDAL  SUBGRAPHS 
AND 

DERIVED  MATRIX  SPLITTINGS 


P.  M.  Dearing 
D.  R.  Shier 
D.  f).  Warner* 

C 1 emson  University 
Clernson,  South  Carolina  29631 


A  linear  system  of  equations  has  a  perfect  elimination  ordering  if 
and  only  if  the  correspond ng  adjacency  graph  is  chordal.  For  systems 
which  do  not  possess  a  perfect  elimination  ordering  it  is  customary  to 
try  to  find  an  ordering  which  minimizes  the  fill-in.  One  such  approach 
is  to  embed  the  corresponding  adjacency  graph  in  a  minimal  chordal 
supergraph. 

Another  approach  for  solving  the  linear  system  is  to  consider  an 
iterative  method  based  on  a  splitting  which  is  defined  by  a  maximal 
chordal  subgraph.  In  this  paper  we  present  an  algorithm  for  finding  a 
maximal  chordal  subgraph  and  explore  the  effectiveness  of  the  derived 
iterative  scheme. 


Session  3A.  12:30  -  12:50  p.m. ,  Monday 


SOME  EXPERIENCES  Of  SOLVINo  '  ARuf  SCALt 
GENERALIZED  EIGENVALUE  PROBLEMS  I ANC/OS  MI THOD 


.Jid'i-<in  Jen 

Computer  Center  of  Academia  Si  ni  r.a 
People's  Republic  of  China 


(Only  title  was  submitted.) 
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A  STATIONARY  ITERATIVE  METHOD  THAT  WORKS  EVEN  FOR 
HERNITIAN  Ax  =  b 


John  de  Pi  II i s* 
University  of  California 
Riverside,  California  92521 


Wi  1  he  1  in  Ni et hammer 
Universitat  Karlsruhe 
D - 7 500  Karlsruhe,  GERMANY 


We  develop  a  theory  for  finding  scalars  ai,  a2,  0:3,  a 4  so  as  to 

produce  the  sequence  defined  by 

yk+2  =  ( aiB+a2)yk+1  +  (a3B+a4)  +  ( aj+a3) A^b,  k  =  0,1,2,... 

where  y0,yi  are  arbitrary.  A  =  A0(I-B)  is  an  m  x  in  non-singular  matrix 
and  matrix  A0  is  easy  to  invert.  The  aq  's  are  chosen  so  that  yk-*x,  and 
Ax  =  b.  Heretofore,  the  theory  of  (stationary)  one  and  two-part 
sequences,  e.g.,  Jacobi,  SO R  methods,  Chebyschev  semi -iterative  method 
(which  is  asymptotically  stationary),  required  o(B),  the  spectrum  of  B, 
to  lie  wholly  on  one  side  or  the  other  of  a  line  passing  through  the 

point  2=1.  For  example,  if  A  is  positive  definite,  A0  =  I, 

a(B)  c=  (-«>,1).  But  if  A  is  hermitian  and  A0  =  I,  o(B)  c=  (1-t,  1+t), 
and  hence  straddles  the  point  z  =  1. 

As  a  special  case  of  our  theory,  we  have  the  theorem: 

Let  A  =  A*  with  spectrum  a(A)  c=  (-t,t).  Then  for  any  y0,yj,  the 

sequence 

yk*2  -  [fl(-yk,i  *  V2>  *  b/^/,R  *  Vi 

always  converges  to  x,  where  Ax  =  b.  The  asymptotic  rate  of  conver¬ 
gence  is  R/\  =  -log(r)  where  r2  =  j  +  (1-K_r)  "72  and  K  is  the  condi¬ 
tion  number  of  A. 
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SPIKE  SELECTION  FOR  LARGE  SPARSE  SETS  OF  NONLINEAR  EQUATIONS 


Arne  Drud 
World  Bank 

Washington,  D.C.  20433 


We  consider  square  nonlinear  systems  with  the  following 
properties : 

a.  The  largest  simultaneous  Block  is  large  (500-3000 
equations) . 

b.  The  system  is  sparse  (3-5  nonzero  derivatives  per 
equation  on  average). 

c.  Most  equations  are  analytically  invertible  w.r.t.  most 
of  its  variables. 

Each  simultaneous  block  of  such  a  system  is  often  solved  with  a 
spike  selection/partitioning/reduction  procedure:  Select  a  set  of 

spike  variables  and  spike  equations  such  that  the  remaining  system  is 
recursive  and  analytically  solvable  for  fixed  spike  variables.  Use  a 
"small  scale"  algorithm  to  solve  the  reduced  system,  i.e.,  to  solve  the 
spike  equations  w.r.t.  the  spike  variables,  implicitly  substituting  out 
all  recursive  variables. 

The  criteria  for  spike  selection  are  to  minimize  the  size  of  the 
reduced  system  and  at  the  same  time  keep  it  wel 1 -condi ti oned. 

The  paper  will  compare  two  spike  selection  algorithms,  one  based 
on  Hel 1 erman-Rari ck  with  post-processing  to  eliminate  noni nverti d1 e 
relationships  in  the  recursive  system  (ref.  swapping  of  triangular 
columns  in  linear  systems),  and  one  that  considers  the  noni n vert i hi  1 i ty 
di rectly . 

The  conditioning  of  the  reduced  system  can,  at  the  cost  of  size, 
be  improved  by  avoiding  certain  unstable  inversions  in  the  recursive 
part.  The  paper  will  report  how  different  i overt i bi 1 i ty  crit-'ria 
influence  the  reduced  size  and  speed  of  convergence  for  some  l.-tr.e 
econonmic  planning  models. 
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SPARSE  MATRIX  TEST  PROBLEMS 


Iain  S.  Duff 
AERE  'Harwell 
Oxon,  0X11  ORA,  ENGLAND 

Roger  G.  Grimes 
John  G.  Lewis 
W.  G.  Poole,  Jr. 

Boeing  Computer  Service,  Inc. 
Tukwila,  Washington  93138 


In  the  spring  of  1932  we  appealed  in  the  IMANA  and  SIGN'JM  news¬ 
letters  for  further  test  matrices  to  augment  collections  at  Harwell  and 
BCS. 

In  this  poster  session,  we  report  on  the  results  of  these  appeals. 
Additionally,  we  indicate  how  we  have  organized  the  current  set  of  test 
matrices  and  how  we  have  classified  the  test  matrices  in  the  collection. 
We  have  designed  our  data  base  for  ease  of  distribution,  to  facilitate 
the  addition  of  further  examples,  and  to  allow  the  same  matrix  to  be 
included  in  more  than  one  cl  ass i f i cat i on . 

We  also  explain  the  mechanism  for  submitting  test  examples  for 
inclusion  in  the  collection  and  for  obtaining  copies  of  the  test 
matri ces . 


Poster  Session.  Monday  and  Tuesday 


THE  MULT  I  FRONTAL  SOLUTION  OF 
UNSYMMETRIC  SETS  OF  LINEAR  EQUATIONS 


Iain  S.  Du  ft 
John  K.  tieid* 

AERL  Harwell 

Oxon,  OX i  1  OKA,  ENGLAND 


We  show  that  general  sparse  sets  of  linear  equations  whose  pattern 
is  symmetric  (or  nearly  so)  can  be  solved  efficiently  by  a  mul tif rontal 
technique.  The  main  advantages  are  tnat  the  analysis  time  is  small 
compared  to  the  f actori zat i on  time  and  analysis  can  be  performed  in  a 
predictable  amount  of  storage.  Additionally,  there  is  scope  for  extra 
performance  during  factori zation  and  solution  on  a  vector  or  parallel 
machine.  We  show  performance  figures  for  examples  run  on  the  IBM  3033 
and  CRAY-1  computers. 


Session  9A.  7:30  -  7:50  p.m. ,  Tuesday 


A  PARTITIONING  APPROACH  10  PUOCLSSINi.  NIIWORKS 
AND  THE  SOLUTION  OF  SPARSI  SYS i I  MS 


Michael  •  ^ ' 

University  o-  as 
Austin,  Texas  /-'I/ 


Ordinary  networks  anise  in  connection  -c'H  r  ne  t  ransportat ion  of 
commodities.  Processing  networss  extend  or. li nary  networks  by  allowing 
proportional ity  constraints  among  certain  Flows. 

When  the  primal  simplex  algorithm  is  applied  to  a  minimum  cost  flow 
problem  for  processing  networks,  the  basis  can  be  partitioned  into  a  set 
of  arcs  which  form  a  so-called  representative  spanning  tree  and  a  set  of 
non-tree  arcs.  This  partitioning  yields  a  working  basis  which  has  a 
lower  dimension  than  that  which  would  result  if  a  (non-representati ve) 
spanning  tree  were  used.  An  implementation  of  the  primal  simplex 
algorithm  has  been  developed  which  maintains  the  basis  in  this  parti¬ 
tioned  form.  Computational  results  will  be  presented  which  show  that 
this  implementation  compares  very  favorably  with  MINOS  on  a  certain 
class  of  processing  networks. 

Sparse  systems  of  equations  can  also  be  modelled  as  processing  net¬ 
works.  It  will  be  shown  that,  in  some  cases,  solving  the  system  involv¬ 
ing  the  working  basis  is  easier  than  solving  the  original  system. 
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MATRIX  MODIFICATION  AND  PARTITIONING 


Albert  M.  •>!>  sa r 
Boeing  Computer  Services  Co  ''panv 
Tukwila,  Washington  T-l:-;-' 


The  contingency  analysis  problem  invo'ves  so! vinj  a  sequence  >f 
systems  of  sparse  equations  where  the  matrices  differ  by  a  .natn>  o* 
low  rank.  These  systems  can  be  solved  by  partitioning  the  matrix  ana, 
confining  the  changes  to  the  border,  or  by  using  the  matrix 
modification  formula  to  obtain  the  so1  .it  ion  of  one  problem  in  terms  or 
the  previous  one. 

Often  the  low  rank  change  has  a  special  form:  the  matrix  can  be 
reordered  so  that  all  changes  are  confined  to  the  lower  right  corner 
block.  If  this  reordering  is  actually  done,  very  little  work  is 
required  in  obtaining  the  solution  of  the  modified  problem  based  on  the 
solution  of  the  original  problen.  Unfortunately,  this  places  a  severe 
constraint  on  the  partitioned  form  which  is  greatly  amplified  when  a 
sequence  of  such  changes  in  different  parts  of  the  matrix  are  desired. 
We  will  show  a  way  of  solving  this  problem,  taking  advantage  of  the 
special  form  of  the  low  rank  change,  but  not  constraining  the  ordering 
of  the  matrix. 
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A  STRUCTURALLY  STABLE  MODIFICATION  OF  I  HE  HELLERMAN-RARICK 
ALGORITHM  FOR  REORDERING  UNSYMMETRIC  SPARSE  MATRICES 


Albert  M.  Erisinan 
Roger  G.  Grimes 
John  G.  Lewis* 

William  G.  Poole,  Jr, 

Boeing  Computer  Services  conran v 
Tukwi la,  Washington  eJI-B 


The  P“  algorithm  of  Hellerman  ,r<-1  -  ar  •.  ,>  -.*s  use::  to  reorder 

unsynimetric  sparse  matrices  in  orb-'  to  r ,  ration  and 

storage  costs  when  solving  sparse  syst--  •,  li"ear  i  ;uat ions. 

Unfortunately,  it  is  not  infallible.  it  •  s  •  f  r  the  algorithm 
can  generate  intermediate  matrices  wrier  -,r~  -,r ;;  ,  ;ul  ar  and 
which  may  lead  to  a  breakdown  in  the  e!  i  r- r. a- 1  .v  :s>. 

In  this  paper  we  present  the  al‘jo<-:tn  m  i  r  for 

ana  explain  several  of  the  problems  which  ,o .  ....r.  i 

several  past  attempts  which  were  unsuccesst  .  : "  i.orre-_:  ’ r  j  tre 

problems . 

We  then  define  a  new  modification  tor  nan.il  •,  r..;  tne  d:  ‘  f  :c..  1  ’  es . 
This  revised  version  of  the  algorithm  will  never  .-rodare  str  ,rj :  ;  . 
singular  intermediate  matrices  if  tne  an  cm  matrix  i  s  ri,* 
structurally  singular.  Test  results  comparing  this  version  with  ot.ners 
will  be  given. 
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SOLVING  SPARSE  STAIRCASE  SYSTEMS 
BY  GAUSSIAN  ELIMINATION 


RoDert  Foarer 
Northwestern  University 
Evanston,  Illinois  60/01 


A  system  of  linear  equations  has  a  "’sparse  staircase"  structure  if 
its  variables  fall  into  a  certain  sequence  of  stages,  such  that  any 
equation  involves  relatively  few  variab1es  from  at  most  two  adjacent 
stages.  Sparse  staircase  systems  commonly  arise  in  solving  large 
multi-period  linear  programs  by  the  simplex  method. 

This  paper  describes  two  approaches  to  solving  sparse  staircase 
systems  by  Gaussian  elimination.  Each  approach  combines  techniques  for 
stable  period-by-period  elimination  of  arbitrary  staircase  matrices 
with  techniques  for  efficient  elimination  of  arbitrary  sparse  matrices. 
The  two  approaches  differ  in  their  choice  of  sparse-el i mi nati on  techni¬ 
ques:  one  uses  "merit"  techniques  that  select  an  attractive  pivot  at 
each  elimination  step,  while  the  other  employs  "bump-and-spi Ke"  techni¬ 
ques  that  permute  the  staircase  matrix  to  an  attractive  form  in  advance 
of  elimination. 

Sparse  staircase  elimination  methods  admit  especially  efficient 
implementations  that  operate  on  only  a  period  or  two  of  the  staircase 
at  a  time.  Their  peri od-by-pen od  elimination  ordering  can  also  be 
advantageous  in  handling  certain  sparse  right-hand-side  and  solution 
vectors,  and  in  computing  solutions  for  subsets  of  periods.  They  may 
also  produce  an  especially  sparse  f acton zat i on  of  the  staircase 
matrix,  although  in  general  tnoy  do  not  find  a  sparser  factori<:ation 
than  general  sparse-el imi nat i on  methods. 
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AN  INCOMPLETE  FACTORISATION  METHOD  FOR 
SINGULAR  IRREDUCIBLE  M- MATRICES 


R.  F .  Funder  1 i c 

Union  Carcide  Corporation,  Nuclear  Division 
Oar  Ridge,  Tennessee  37030 

R.  J.  ^lemmons 

North  Carolina  State  University 
Raleigh,  Nortn  Carolina  27650 


Markov  queueing  networks  can  give  rise  to  very  large,  sparse, 
irreducible,  singular  (zero  column  sums)  M-natnces  A.  The  stationary 
probability  distribution  vector  is  the  solution  to  a  homogeneous  system 
of  linear  equations  Ax  =  0.  Certain  economic  petroleum  reservoir  and 
discrete  Neumann  problems  give  rise  to  related  systems.  The  splitting 
A  =  M  -  N  with  a  matrix  M  having  symmetric  zero  structure  can  be  shown 
to  be  a  regular  splitting.  A  sparse  LU  factori zati on  of  a  symmetric 
permutation,  PMPT,  is  obtained  using  a  standard  symmetric  ordering  such 
as  minimum  degree.  No  pivoting  for  stability  is  necessary.  Splitting 
strategies  including  those  that  have  the  larger  1  a  j  j  j  in  M  will  be 
discussed. 
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ROW  ORDERING  SCHEMES  FOR  SPARSE  GIVENS  TRANSFORMATIONS 


-Man  George 

University  of  Waterloo 
Waterloo,  Ontario,  Canaaa  N21  3G1 

Joseph  Liu* 

York  Unveristy 

Downsview,  Ontario,  Canaaa  M3J  1 P 3 


Let  A  be  an  m  by  n  matrix,  in  >  n,  and  let  P;  and  P-  be  permutation 
matrices  of  order  m  and  n  respecFi  vely ,  Suppose  p  AP  is  reduced  to 

Rs 

upper  trapezoidal  form  ('  '  using  Givens  rotations.  The  sparsity 

structure  of  R  depends  only  on  p2.  For  a  given  P2,  the  number  jf 

arithmetic  operations  required  to  compute  R  depends  on  P,.  In  this 
paper  we  consider  row  ordering  strategies  that  are  appropriate  when  P, 
is  obtained  from  nested  dissection  type  orderings  of  A'A.  Recently  it 
was  shown  that  the  so-called  "width-21  nested  dissection  orderings  of 
aTa  could  be  used  to  simultaneously  obtain  good  row  and  column  order¬ 
ings  for  A.  In  this  paper  we  show  that  the  conventional  (width-1) 

nested  dissection  orderings  can  also  be  used  to  induce  good  row 

orderings.  Our  analysis  employs  a  bi-partite  graph  model  of  Givens 
rotations  applied  to  A,  similar  in  some  respects  to  the  bi-partite 
graph  models  of  Gaussian  elimination  developed  by  Gilbert  and  Mauli^o. 


Session  88.  12:10  -  12:30  p.m. ,  Tuesday 


SOLUTION  OF  SPARSE  UNDERUETERMI NEU 
SYSTEMS  OF  LINEAR  EQUATIONS 


Alan  Ueorye 
Esmond  Ny* 

University  of  Water'00 
Waterloo,  Ontario,  Canada  N21  3'ii 


In  this  paper  we  consider  the  problem  of  computing  tne  minimal 
^-solution  to  a  consistent  underdetermined  linear  system  Ax  -  b,  where 
A  is  in  by  n  with  m  <  n  .  The  method  of  solution  is  to  reduce  A  to 
lower  trapezoidal  form  [L  0]  using  orthogonal  transformations ,  where  L 
is  m  by  m  and  lower  triangular.  The  method  can  be  implemented  effi¬ 
ciently  if  the  matrix  AA^  is  sparse.  However,  if  A  contains  some  dense 
columns,  AA^  may  be  unacceptably  dense.  We  will  present  a  method  for 
handling  these  dense  columns.  The  problem  of  solving  a  rank-def i ci ent 
under-determined  system  will  also  be  considered. 
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NESTED  DISSECTION  WITH  PARTIAL  PIVOTING 


John  R.  Li i  Inert* 
Cornel  1  Un 1  vers i ty 
Ithaca,  New  Yorx  14853 

Robert  Schreiber 
Stanford  University 
Stanford  California  94305 


When  a  sparse  system  of  linear  equations  is  solved  by  Gaussian 
elimination,  many  of  the  zero  coefficients  often  become  nonzero.  One 
strategy  for  limiting  this  fill-in  is  nested  dissection,  which  is 
applicable  to  many  two-dimensional  finite  element  and  finite  difference 
problems.  Alan  George  invented  this  algorithm,  and  Lipton  and  Tarjan 
extended  it  to  give  fill-in  that  is  within  a  constant  factor  of  minimum 
for  any  system  whose  coefficient  matrix  represents  a  planar  or 
near-planar  graph. 

Most  analyses  of  nested  dissection  avoid  numerical  problems  by 
assuming  that  the  matrix  of  coefficients  is  symmetric  and  positive 
definite.  What  happens  if  the  matrix  is  not  so  well  behaved?  This 
talk  will  present  an  algorithm  that  applies  nested  dissection  while 
performing  partial  pivoting  to  control  numerical  stability.  We  shall 
analyze  this  "dissection  pivoting"  algorithm  by  using  a  bipartite  graph 
model.  For  a  large  class  of  matrices,  including  those  that  represent 
planar  graphs  of  bounded  degree,  the  fill-in  from  this  algorithm  is 
within  a  constant  factor  of  minimum.  The  constant  is  larg_,  and  we 
shall  discuss  what  can  be  done  to  make  the  algorithm  practical. 


Session  9A.  7:50  -  8:10  p.m. ,  Tuesday 


IMPLEMENTATION  OF  A  LARGE-SCALE  LINEAR 
PROGRAMMING  PACKAGE  ON  A  MINICOMPUTER 


James  E.  Giles 
Tennessee  Valley  Authority 
Norris,  Tennessee  378?8 


Implementation  of  the  linear  programming  portion  of  the  software- 
package  MINOS  on  the  KP1000F  minicomputer  is  discussed.  Memory 
restrictions  are  circumvented  using  the  extended  memory  addressing 
capabilities  of  the  HP1000F.  An  application  of  the  implemented  package 
to  a  water  resource  management  problem  is  described. 
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ANALYSIS  OF  A  MULTIGRID  METHOD  AS  AN  ITERATIVE  TECHNIQUE 
FOR  SOLVING  LINEAR  SYSTEMS 


Anne  Greenbaum 

Lawrence  Livermore  National  Laboratory 
Livermore,  California  945S0 


A  general  class  of  iterative  methods  is  introduced  tor  solving 
symmetric,  positive  definite  linear  systems.  These  methods  use  two 

different  approximati ons  to  the  inverse  of  the  matrix  of  the  problem, 
one  of  which  involves  the  inverse  of  a  smaller  matrix.  It  is  shown 
that  the  methods  of  this  class  reduce  the  error  by  a  constant  factun 
at  each  step  and  that  under  "ideal"  circumstances  this  constant  is 
equal  to  ( -l)/( «'  +1 ) ,  where  is  the  ratio  of  the  largest  eigen¬ 

value  to  the  (J+l)st  eigenvalue  of  the  matrix,  J  being  the  dimension 
of  the  smaller  matrix  involved.  A  multigrid  method  is  presented  as 

an  example  of  a  method  of  this  class,  and  it  is  shown  that  while  the 
multigrid  method  does  not  quite  achieve  this  optimal  rate  of  conver¬ 
gence,  it  does  reduce  the  error  at  each  step  by  a  constant  factor 

independent  of  the  mesh  spacing  h. 
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PROGRAM  FOR  THE  COMPARISON  OF  REORDERING  ALGORITHMS  FOR  THE 
SOLUTION  OF  UNSYMMETRIC  SPARSE  SYSTEMS  OF  EQUATIONS 


Roger  G.  Grimes 
Jonn  G.  Lewis 
Wi 1 1 i am  G.  Pool e,  Jr. 

Boeing  Computer  Services  Company 
Tukwi la,  Washington  98188 


The  authors  have  developed  a  program  to  monitor  the  effectiveness 
of  various  reordering  algorithms  for  the  solution  of  unsymmetric 
systems  of  equations  over  a  broad  collection  of  test  problems.  The 
objective  is  to  select  implementation-independent  measures  which  can 
lead  to  the  classification  of  matrices  into  categories  where  a 
particular  reordering  is  most  effective.  Measures  chosen  include  the 
amount  of  fill  generated,  operation  counts,  and  numerical  stability  in 
the  factorization  of  the  permuted  matrix. 

A  'arge  collection  of  test  matrices  from  diverse  applications  is 
being  assembled  to  be  used  in  analyzing  matrix  patterns.  A  description 
of  the  design  and  implementation  of  the  program  along  with  samples  of 
preliminary  results  will  be  presented. 
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A  NEW  VERSION  OF  TARJAN’S 
STRONG  CONNECT  ALGORITHM 


Fred  'j.  uustavsor 
IBM  T.  J.  Watson  Researcn  Center 
Yorktown  ^eights,  New  Yor*  10593 


We  describe  a  new  version  ct  ar  ;  algorithm  for  finding  the 
strong  components  of  a  directed  graph  n  wtn  n  nodes  and  N  edges.  The 
new  version,  called  8LTF  (Bloc*  Lower  'n a"  lar  Form),  is  presented  in 
a  structured  high  level  ton  tnat  is  trans'sfed  into  an  efficient 
FORTRAN  code.  This  version  runs  faster  tha°  a  previous  version  called 
STCO.  This  algorithm  is  closely  related  to  the  algorithm  MC13D  of  Duff 
and  Reid,  which  is  also  an  implementation  of  STRONG  CONNECT.  Algorithm 
BLTF  is  superior  to  MC13D  in  the  following  ways:  it  executes  faster  and 
uses  less  storage;  it  computes  exactly  the  same  output  as  Tarjan's 
algorithm;  its  translation  into  FORTRAN  is  completely  structured;  its 
inner  loop  is  more  efficient  in  that  it  asxs  half  as  many  questions  as 
MC130  as  N  -  n2. 
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AN  ALGORITHM  FOR  COMMUTING  A  FULL  ASSIGNMINT 
FOR  A  SPARSE  MATRIX 


Fred  G.  Gustavson 
IBM  T.  J.  Watson  Research  Center 
Yorktown  Heights,  New  York  105R8 


We  describe  an  efficient  implementation  of  the  Assign  Row 
Algorithm  of  Gustavson.  This  algorithm  finds  a  maximal  assignment  for 
an  arbitrary  sparse  (0-1)  matrix  and  is  based  on  the  algorithm  of 
M.  Hall.  This  algorithm  is  closely  related  to  algorithm  MC21A  of  Duff, 
which  is  also  an  implementation  of  Hall's  algorithm.  Two  versions  of 
Assign  Row  are  discussed.  The  first  computes  all  cheap  assignments 
before  starting  any  depth  first  searches.  The  second  implements  Assign 
Row.  Both  algorithms  are  completely  structured  and  efficiently 
translated  into  FORTRAN.  The  first  algorithm  almost  always  computes 
faster  than  the  second  and  both  algorithms  almost  always  computes 
faster  than  MC21A. 
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AN  IMPROVED  ASSIGNMENT  ALGORITHM 


Fred  G.  Gustavson* 

Gary  D.  Hachtel 

IBM  T.  J.  Watson  Research  Center 
Yorktown  Heights,  New  York  10598 

Aram  K.  Kevorkian 
General  Atomic  Company 
San  Diego,  California  92138 


The  Assign  Row  Algorithm  of  Gustavson  uses  the  depth-first  search 
technique  to  stretch  an  assignment  in  an  undirected  bipartite  graph. 
An  efficient  implementation  of  this  algorithm  is  the  fastest  method  for 
finding  an  assignment  for  almost  all  graphs. 

We  present  a  new  and  improved  version  of  the  Assign  Row  Algorithm 
which  performs  equally  well.  However,  on  the  small  subset  of  all 
graphs  where  the  Assign  Row  Algorithm  performs  poorly,  the  new 
algorithm  runs  up  to  an  order  of  magnitude  faster.  The  main  feature  of 
the  new  algorithm  is  the  introduction  of  the  efficient  methods  in  the 
depth-first  search  procedure  so  that  negative  results  in  earlier  depth- 
first  searches  can  be  used  to  reduce  present  and  future  depth-first 
searches.  The  new  algorithm  sheds  insight  into  the  problem  of  finding 
the  block  lower  triangular  form  of  a  sparse  matrix  as  follows.  A 
partial  assignment  gives  rise  to  a  'partial  directed'  graph  and  depth- 
first  search  can  be  used  to  quickly  explore  this  graph.  Failure  to 
stretch  an  assignment  in  this  graph  leads  to  partial  information  on  the 
block  lower  triangular  form.  Hence  these  parts  of  the  graph  need  not 
be  explored  in  later  depth-first  searches.  Furthermore,  the  basic 
technique  used  in  the  two  stage  procedure  of  finding  an  assignment 
followed  by  the  strong  components  determination  of  the  associated 
directed  graph  merges  into  one  idea  of  repeated  depth-first  search  on 
a  partial  directed  graph. 
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A  1-PASS  PROCEDURE  FUR  MAXIMAL  ASj IGNMFN T 
AND  FINDING  THL  BLT  FORM 


odfy  ;.i .  HaOite  i  * 

Fred  u.  uustavsoo 
IBM  r.  J.  Watson  Research  Center 
for«,town  Heights,  New  v')r,  Cjv<  > 


Aran  n.  Kevorkian 
General  Atomic  Company 
San  Diego,  California  r’IS-; 


A  procedure  is  proposed  for  obtaining  noth  a  maxima'  assignment 
and  the  canonical  blocx  triangular  form  of  a  sparse  matrix.  In 
contrast  to  the  conventional  2-pass  approach,  the  proposed  unified 
procedure  incorporates  block  identification  tests  into  the  assignment 
algorithm.  Thus  blocks  of  the  BLT  form  are  discovered  while  in  the 
process  of  solving  the  0'|Vj:.?  (|V|  +  jE|)  •  maximal  assignment  problem. 

This  procedure  is  shown  to  simplify  and  decompose  the  maximal  assign¬ 
ment  problem  by  virtue  of  removing  large  row  and  columns  blocks  from 

further  consideration.  Each  irreducible  blocx  is  identified  once  only. 
The  accumulated  cost  of  the  test  is  linear  n.e.,  (J(  j  V  I  +  |  E  j )  ).  The  new 
procedure  is  based  upon  OF  S  of  the  corresponding  bipartite  graph. 

Alternating  edges  in  the  palm  tree  thus  obtained  give  a  partial 
assignment  of  greater  cardinality  than  that  obtained  by  previously 
published  "cheap  assign"  methods.  Since  the  running  time  of  maximal 
assignment  algorithms  depends  chiefly  on  the  number  of  subsequent 
"stretching"  DFS  passes,  this  reduces  the  expense  of  the  overall 

procedure. 

The  partial  assignment  makes  the  hi  part  1 te  graph  partially 

directed.  When  backtracking  during  PFS  of  the  bipartite  graph,  we 

perform  tests  analagous  to  Tarjan's  transitive  tests  for  bi connected 
components  of  the  bipartite  graph  and  strong  components  of  the  directed 
subgraph  (both  row-wise  and  col  -wise) .  Tnese  teds  can  identify 

irreducible  blocks  o*  the  BLT  tor"‘  'or  set-,  of  ...ch  b  1  uo*  s  '•  p  •'  1  o  r  to 

completing  the  ass  i  unmeet .  we  )’,n  -  1™  >.»rf  am  ;,m-  •’mi  mum 

cover"  tests  which  can  :  h»*»t  i  fy  So  *•  o’  e  ■  s  the  ,*  W  •  r.Jt,<.>«v’  of 

backtracking.  rhese  test  >  re  i  i  ’  re  less  t  *t  i  ;ml  mpense,  smee 

they  are  not  transitive. 

We  present  some  stat  '  d  1  c a  1  d  it  i  whi  *>  *  ’•  *’v  i  rm -m  , 

might  be  favorable  when  up!  mi  t  ;  ’  .  •  : 1  ..-■■■  ■  :  -  *  r 

reduci bi ! i ty  . 
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A  SPARSE  LINEAR  PROGRAMMING  SUBPROGRAM 


S.  J.  Hanson 
K.  L.  Hiebert* 

Sandia  National  Laboratories 
Albuquerque,  New  Mexico  87135 


We  will  describe  a  subprogram,  SPLP,  for  solving  linear 
programming  problems.  The  package  of  subprogram  units  comprising  SPLP 
is  written  in  FORTRAN  77. 

The  subprogram  SPLP  is  intended  for  problems  involving  at  most  a 
few  thousand  constraints  and  variables.  The  subprograms  are  written  to 
take  advantage  of  sparsity  in  the  constraint  matrix. 

A  very  general  problem  statement  is  accepted  by  SPLP.  It  allows 
upper,  lower,  or  no  bounds  on  both  the  variables  and  the  constraints. 
Both  the  primal  and  dual  solutions  are  returned  as  output  parameters. 
The  package  has  many  optional  features.  Among  them  is  the  ability  to 
save  partial  results  and  then  use  them  to  continue  the  computation  at  a 
later  time. 
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SPARSE  ORTHOGONAL  SCHEMES  TOR  SI  RUG  T  URAL  OP !  IMITATION 
USING  Till  LORO.  METHOD 


M 

Union  Carbide  corp_. 

Od  k  R  ’  d  q  e , 
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R.  C.  ward 

Union  CarDide  Corporation,  Nuclear  ••  .■  :  s  i  on 
Gas  Rid ye,  Tennessee  3. Rio 


The  force  method  enab'es  one  to  calculate  the  system  force  vector 
F  acting  on  a  structure,  and  tne  associated  stresses  and  strains, 
without  reassembling  and  refactoring  the  overall  stiffness  matrix  at 
each  step  of  the  optimization  process.  The  me  mod  consists  of  two 
parts:  (1)  computation  of  a  matrix  B  whose  col  mins  form  a  basis  of 

the  null-space  of  the  equilibrium  matrix  E,  which  is  m  *  n  of  rank  m, 
along  with  a  particular  solution  S  to  the  underdetermined  system 
ES  =  P,  where  P  is  the  set  of  external  loads  on  tne  structure,  and 

{?.)  tne  solution  of  the  least  squares  problem  nun  if  (3x  +  S)  r?,  where 

y 

f  -  f 1  f  is  the  element  flexibility  matrix  which  changes  at  each 
optimization  step.  In  this  case  f  -  Bx  +•  S  satisfies  the  principal  of 
complementary  potential  energy  -  that  E ■ f -  is  r’nrnal  over  all  F  such 
that  EF  -  p. 

I  n  this  talk,  an  orthogonal  iecomposi  t  i  c-  ice  c  for  solving  part 
(1),  called  the  turnbacx-QR  method,  *nic.h  .-reserves  me  banded  struc¬ 
ture  of  E  in  computing  B,  is  presented.  Some  ;ir«* I  imirary  comparisons 
are  made  with  the  f  urnba; «. -Eli  net  ho;:  oa-m-i  uoon  e!  ’mi nat i on,  r  ver’ 
recently  by  Kaneko,  et  a  1  ’leBBl.  -'-iv.  s  .  •••r.-nts  .r.'  •  on  fv 
ise  of  the  1  i  near  ly  c  •••• :  s  ;  ,  y  ,•  •■•••  no-.'  e  t  or 

computing  x  ’  c  >7)  at  mm  o  i  .»  r.  •  . 


-  \?  :bii  fi.m.  ,  Monday 
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INCOMPLETE  METHODS  FOR  SOLVING  ATAx  -  b 
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Bellas*.  Northern 
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Much  interest  has  recently  been  shown  in  incomplete  f actor i cat i on 
methods  for  solving  sparse  linear  simultaneous  equations.  They  hdve 
the  advantage  of  avoiding  most  or  all  of  the  fill-in  associated  with 
pure  elimination  methods  and,  although  iterative  in  nature,  tney 
exhibit  much  better  convergence  rates  man  do  classical  iterative 
methods . 

In  this  paper  incomplete  methods  are  considered  tor  tne  solution 
of  A^Ax  =  b  where  A  is  sparse.  One  possibility  is  to  construct  the 
product  of  AtA  and  then  use  an  ICCG  algorithm  for  tneir  solution. 
Alternatively  it  is  possible  to  perform  incomplete  orthogonal  decompo¬ 
sitions  of  A  by  modifying  either  the  Gram-Schmi dt  or  the  Givens  rota¬ 
tion  methods.  These  decompositions  will  yield  incomplete  trianguiar 
factors  for  A^A  which  differ  from  tne  one  obtained  by  ICCG,  but.  whicn 
give  a  similar  iterative  phase  for  the  solution  process. 

Implementations  of  each  of  these  three  methods  are  briefly 
described  and  some  numerical  cornua r  i  so-s  cade  using  examples  from 
structural  analysis.  The  ICCG  method  has  *v  advantage  that  it  is  easy 
to  implement  using  passing  schemes  me  elements  of  the  -.parse 

matrices.  The  orthogonal  deco: pos’t  ion  . is .  particularly  the  one 

using  Givens  rotations,  do  atueir  to  n  a  s  ,i  •«whut  netter  numer  ica  i 
efficiencies  when  counting  only  non-zero  arithmetical  operations,  but 
it  is  not  obvious  how  them  os>f  i,  ••  ,  u.jses  M  he  el *■  i  r  f  ent 

implemented  in  a  sparse  store. 
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A  PRODUCTION  EIGENANALYSIS  SYSTEM  FOR 
LARGE,  GENERALIZED  SYMMETRIC  PROBLEMS 


Paul  S.  Jensen 

Lockheed  Palo  Alto  Research  Laboratory 
Palo  Alto,  California  94304 


A  system  of  programs  called  BES  (Basic  ti genana iysi s  System)  for 
determining  a  few  eigenpairs  of  problems  of  the  form 

Ax  =  Mx  \  ( 1 ) 

has  been  developed  for  production  engineering  analysis.  A  and  M  are 
typically  large,  real,  symmetric  and  sparse  matrices  for  which  there 
exist  scalars  a  and  u  such  that  a A  -  uM  is  positive  definite.  The 
eigenvalues  of  such  problems  are  real  and  the  desired  eigenpairs  are 
typically  specified  by  giving  the  limits  of  a  desired  spectral  range 
(or  section). 

The  system  is  organized  as  independent  executable  programs  that 
operate  in  a  chained  (pipeline)  fashion,  communicating  via  a  database 
system.  For  each  problem  analyzed,  the  chaining  details  are  developed 
dynamically  by  the  system  to  suit  the  particular  problem  and  the 
results  are  deposited  in  a  global  database. 

For  the  solution,  a  series  of  problems  of  the  form 

(A  -  oM)y  -  (aA  -  ufl)y  y  (2) 

is  solved,  where  a  and  g  are  selected  once  for  positive  definiteness 
and  "shift  point"  a  varies.  Each  problem  in  form  (2)  is  transformed  to 
a  classical  form  and  solved  using  a  block  Lanczos  algorithm,  orginiated 
by  Parlett  and  Scott,  that  is  specifically  designed  for  spectral 
sections.  The  transf ormations  require  one  factorization  of  xA  -  uM  and 
one  factori zation  of  A  -  M  for  each  value  of  a. 

Each  shift  point  a  is  at  the  center  of  a  subsection  (interval) 
within  the  range  specified  by  a  user.  For  each  subsection,  the  Lanczos 
analysis  either  determines  all  of  the  eigenvalues  contained  therein  or 
it  determines  several  in  the  vicinity  of  the  shift  point.  In  t lie 
former  case,  no  further  analysis  of  the  subsection  is  performed.  In 
the  latter  case,  two  new  subsections  (at  the  lower  and  upper  extremes 
of  the  given  subsection)  are  generated  ana  subsequently  analyzed.  The 
usual  checks  for  completeness  of  coverage  are  utilized  and  refinements 
via  Ritz  projections  are  made  when  needed. 

A  variety  of  probl  e  ns  f  -'om  structural  v  i  ■'  r  it ;  >r  and  'nn.x  !  vi  j 
analysis  have  been  solved  using  BEs.  It  igear,  *  r  .m  r  typical 

analyses  for  about  one-sixth  of  tee  cost  aMr,1v,*>-a  to  previ  o..s 
programs  based  on  s  in  > '  taneuos  invor  ■  ■  t  ♦  o  ••  a  *  -a  . 
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RESEMBLING  FHOSt  OF  KAFR  ICES  IN  A  SMAUTK  0  IMF  NS  ION 


on  1 


For  a  2-by-2  component -w i oe  -  j t- :  i  x ,  both  eigenvalues 

are  real  and  the  spectral  radios  is  easi  iv  calc^’dte.l  or  estimated. 
For  an  n-by-n  tri -diagonal  nonnegative  -;a t r 1 x  these  observations  gener¬ 
alize:  it  is  we  1 1  known  that  the  ^i  qenva  I  ..nt  or-?  all  real,  and  it  nay 

be  shown  that  the  specta!  radios  1 s  o,;  the  order  of  'in  fact,  between  1 
and  2  times)  the  largest  of  the  spectral  ra-li  i  of  the  2-by-2  principal 
submatrices  occurring  consecutively  along  the  diaoonal .  These  facts 
support  the,  perhaps  often  made,  meta-observat ion  that  properties  of 
tri -diagonal  matrices  are  very  much  like  those  of  ?-by-2  matrices. 
However,  both  these  facts  generalize  considerably  further  in  addressing 
the  issue  of  what  it  is  about  tne  sparsity  pattern  of  a  large  matrix 
with  many  zeroes  which  make  it  "like"  a  matrix  in  a  much  smaller  dimen¬ 
sion.  If  the  longest  simple  circuit  in  the  usual  directed  graph  of  an 
n-by-n  matrix  is  k  (<n),  we  snail  say  that  it  belongs  to  the  class 

M  ,  .  For  example  M  „  are  the  ni 1  potent  essentially  triangular 

n  y  K  n  ,  ij 

matrices  M  ,  are  the  remaining  essentially  trianguia-  matrices;  and 

n ,  i 

M  -  includes  (but  is  not  restricted  to)  tie  nontrivial  tri -diagonal 
n  ,2 

matrices.  We  present  here  some  recent  results  about  ways  m  which 
^  is  more  like  the  k-by-k  matrices  tnan  tne  n-by-n  matrices,  regard¬ 
less  of  the  size  of  n.  These  include  sr*  :  t r  a  i  properties  c  ron'-vga- 
tive  matrices  and  general  matrices,  H  -  ;tan i  '  '  t v ,  neter  - 1 ranta !  cal.  na¬ 
tions,  and  others,  nr  ten,  attempts  ire  .r.j  ie  to.  general  ize  tri  -di  agora  1 
results  to  penta-di  agonal  matrice-,,  ef.  . ,  ,-,t  as  ?-ny-2  results  might 

be  generalized  to  the  ll-by-3  case.  Huwov  ;r ,  i f  tne  tri -«!  i  agcna1  tact 
is  actually  an  Mt)  tact,  general  izat  ion  u.  y  (the  onnte-di  ago-nai  s 

are  genera!!/  in  M  ,  -:>vy  be  -  ;.-e  ‘  •  ■  t  *  .  •  . 
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AN  EFFICIENT  ALGORITHM  FOR  THE  FACTORIZATION 
OF  POSSIBLY  RANK  DEFICIENT  MATRICES 


M.  S.  Kamel* 

K.  Singhal 

University  of  Waterloo 
Waterloo,  Ontario,  Canada  N2L  301 


Rank  deficient  matrices  arise  due  to  the  problem's  inherent  struc¬ 
ture,  or  due  to  data  uncertainties  or  perturbations  'nduced  by  round¬ 
off  errors.  It  is  desired  to  have  algorithms  whicn  can  detect  the 
possibility  of  rank,  deficiency  of  a  matrix  in  order  to  -eplace  it  by  an 
actual  rank  deficient  matrix.  Lawson  and  Hanson  propose  an  algorithm 
that  uses  Householder  transf ormat i ons  with  column  interchange  to  trian- 
gularize  the  matrix  and  determine  its  pseudorank.  Their  algorithm  is 
mainly  focused  on  solving  least  squares  problems. 

In  this  paper  we  propose  a  similar  algorithm  to  factorize  the 
matrix  using  Householder  transformations  also  with  column  i nterchanyes . 
However,  the  new  algorithm  avoids  operations  on  those  parts  of  the 
matrix  that  will  eventually  be  discarded  due  to  rank  deficiency.  The 
algorithm  offers  computational  savings  in  the  general  case,  but  its 
major  advantage  is  that  it  can  be  applied  to  sparse  matrices  without 
the  risk  of  ruining  the  sparsity  pattern  that  is  usually  experienced 
when  applying  Householder  transf ormat l ons  to  sparse  matrices. 

In  this  paper  we  present  the  algorithm,  operation  count  analysis, 
applications  to  least  squares  problems,  and  other  appl i cati ons . 
Numerical  experiments  and  examples  will  Dt  reported. 


Session  8B.  12:50  -  1:10  p.m.,  Tuesday 


A  THREE  DIMENSIONAL  rlNIU.  ELLMLNI 
MAGNETOTELLURIC  MODELING  PROGRAM 


w.  P.  .vamp 

Phillips  Petroleum  Company 
Bartlesville,  Oklahoma  74004 


A  three-dimensional  modeling  program  for  maynetotel 1 ur 1 cs  mas  been 
written  using  the  finite  element  method.  The  program  has  been 
successfully  implemented  on  the  CRAY  and  IBM  3U33  computers.  We  will 
discuss  the  experience  gained  by  this  exercise  both  froni  a  numerical 
analysis  and  geophysical  point  of  view.  Example  models  together  with 
comparisons  to  known  models  will  be  presented. 

The  use  of  various  sparse  matrix  packages  in  solving  this  problem 
will  be  discussed.  Of  particular  interest  will  be  the  application  of 
complex  versions  of  the  algorithms  on  the  CRAY  computer. 
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THE  TURN-BACK  LU  PROCEDURE  FOR  COMPUTING 
A  SPARSE  AND  BANDED  BASIS  OF  THE  NULL  SPACE 


Ken  KanexO 

University  of  Wisconsin 
Madison,  Wisconsin  b3/"i)6 


This  talk  is  concerned  with  the  problem  of  conput 1  n-.j  ti  basis  of 
the  noil  space  of  a  given  large  sparse  and  handed  matrix  A  with  full 
row  rank.  The  problem  is  not  only  one  of  the  basic  problems  in  linear 
numerical  analysis,  but  also  has  many  applications  of  practical 
importance.  Two  significant  such  applications  are  (i)  structural 
analysis,  and  (ii)  a  new  implementation  of  the  Ellipsoid  Algorithm  for 
solving  a  linear  system,  or  a  pair  of  dual  linear  programs. 

Determining  any  basis  of  the  null  space  of  a  given  matrix  A  is  not 
a  particularly  difficult  task  and  there  are  many  existing  methods  to  do 
the  task.  Existing  methods,  however,  tend  to  destroy  either  sparsity 
and/or  the  banded  nonzero  pattern  of  A.  In  a  recent  paper,  we  have 
shown  that  a  method  which  we  call  the  Turn-Back  LU  Procedure  is 
extremely  effective  in  computing  a  basis  of  the  null  space  of  A,  while 
retaining  the  sparsity  and  bandedness.  We  shall  show  how  this  Turn- 
Back  LU  Procedure  works  and  demonstrate  its  superiority  over  a  few 
existing  methods  by  means  of  several  examples  arising  from  structural 
analysis. 


Session  3B.  12:10  -  12:30  p.m.  ,  Monday 


FAST  DIRECT  METHODS  FOR 
OTHER  SPARSE  MATRIX  PROBLEMS 


Linda  Kaufman 
8e 1 1  Laboratories 
Murray  Hill,  New  Jersey  07974 


rn 

An  n*n  matrix  A  is  separable  if  it  can  be  written  as  A  =  J  B-j  ®  C-j 

i  =1 

where  the  dimension  of  the  B-j  ‘s  is  greater  than  1  and  there  exist 
matrices  Q  and  Z  such  that  for  all  i,  QB-jZ  are  diagonal  matrices.  Fast 

direct  methods  based  on  matrix  separability  have  been  used  to  solve  the 
linear  system  arising  for  a  finite  difference  discretization  of 
Poisson's  equation.  We  will  show  how  matrix  separability  can  be  used 
to  formulate  efficient  algorithms  in  two  other  contexts.  The  first 
system  arises  when  applying  Galerkin's  method  to  solve  separable 
elliptic  p.d.e.'s  on  a  rectangle  while  using  a  tensor  product  B-spline 
basis.  Using  the  separability  of  the  system  produces  a  dramatic 
decrease  in  time  and  space  requirements  over  traditional  sparse  direct 
methods.  The  second  situation  arises  when  determining  the  probability 
distribution  during  a  quequeing  analysis  of  a  network.  Because  many  of 
the  diagonal  submatrices  of  the  matrix  defined  by  the  Kolmogorov 
balance  equations  are  often  separable,  a  fast  block  ite.ation  technique 
can  be  formulated. 
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A  NtW  THEORY  ON  PERMUTING  MATRICES 
10  BLOCK  TRIANGULAR  FORM 


"‘r,)::  f. .  hevor«i  an* 
i^iera'i  At oimc  Company 
-cn  California  92138 


’  rV'l  i .  .r.i s  t  a  v so n 
Cary  D.  Hachtel 

1 8M  T.  Watson  Research  Center 
for* town  Heights,  'Yew  York  10598 


In  this  paper  we  introduce  the  concept  of  block  transversal  in  a 
matrix.  We  use  this  concept  to  establish  a  senes  of  results  on  block 
triangular  matrices  with  fully  irreducible  (indecomposable)  diagonal 
blocks.  Our  results  include  a  new  theorem  on  the  uniqueness  of  the 
block  triangular  form.  Additionally,  we  employ  the  concept  of  block 
transversal  to  generalize  a  generally  accepted  two-step  procedure  for 
computing  the  block  triangular  form  of  a  square  matrix. 


Session  7A.  10:40  -  11:00  a.m. ,  Tuesday 


BLOCK  ITERATIVE  METHODS  WITH  AGGREGATION  FOR  SOLVING 
NEARLY  COMPLETELY  DECOMPOSABLE  MARKOV  CHAINS 


Robert  Koury 
David  F .  McAl 1 i ster 
William  J.  Stewart* 

North  Carolina  State  University 
Pa  1 e i gh ,  North  Carolina  27650 


Iterative  methods  have  long  been  used  to  obtain  the  stationary 
probability  vector  of  Markov  chains.  These  chains  give  rise  to 
stochastic  matrices  which  are  often  very  large  and  extremely  sparse. 
Furthermore,  in  most  practical  applications  the  matrices  possess  a 
distinctive  non-zero  structure.  The  Markov  chain  is  said  to  be  nearly 
completely  decomposable  when  it  is  possible  to  symmetrically  permute 
the  matrix  so  that  the  probability  mass  is  concentrated  into  diagonal 
blocks;  the  non-zero  elements  of  the  off-diagonal  blocks  being 
relatively  small  in  magnitude.  Under  these  circumstances  the  rate  of 
convergence  of  the  usual  iterative  methods  is  so  slow  that  they  are  of 
no  practical  value  and  analysts  have  been  obliged  to  turn  to  decomposi- 
tional  and  aggregative  techniques  to  determine  approximate  solutions. 
The  accuracy  of  the  approximati on  obtained  is  in  general  proportional 
to  the  degree  to  which  the  probability  mass  is  concentrated  into  the 
diagonal  blocks. 

Recent  research  has  lead  to  a  number  of  methods  in  which 
aggregation  is  incorporated  into  iterative  procedures  to  enable  exact 
solutions  to  be  computed  efficiently.  In  this  paper  we  distinguish  two 
different  aggregation  techniques  and  discuss  their  applicability.  We 
show  how  these  techniques  may  be  combined  with  the  group  iteration 
methods  of  Gauss-Seidel  and  successive  overrelaxation,  and  we  present 
the  results  which  were  obtained  from  a  number  of  important  examples. 


Session  3A.  12:10  -  12:30  p.m. ,  Monday 


A  COMPARISON  OF  TWO  SPARSE  MATRIX 
PROCESSING  TECHNIQUES  FOR 
STRUCTURAL  ANALYSIS  APPLICATIONS 


William  A.  Loden 

Lockheed  Missiles  and  Space  Company,  Inc, 
Sunnyvale,  California  94086 


The  performance  of  algorithms  based  on  two  techniques  for  the 
solution  of  large  sparse,  symmetric  linear  equation  systems  of  the 
type  that  arise  frequently  in  structural  analysis  are  studied. 
Representative  test  problems,  for  simple  and  complex  structural 
configurations  treated  with  the  finite  element  method,  are  solved 
with  two  implementations  of  Cholesky  method  profile  algorithms  and 
with  the  SPSYST  sparse-matri x  system  package  developed  by  J.  K.  Reid. 
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CONJUGATE  GRADIENT  IMPLIES  NORMAL 


Tom  Manteuffel* 

Vance  Faber 

Los  Alamos  National  Laburatory 
Los  Alamos,  New  Mexico  87545 


A  conjugate  gradient-like  iteration  for  general  linear  systems  has 
long  been  sought.  The  existence  of  such  a  method  depends  upon  how 
"conjugate  gradient-like"  is  defined.  One  can  always  form  the  normal 
equations  and  apply  the  conjugate  gradient  method  for  symmetric 
positive  definite  systems.  We  consider  gradient  methods  as  defined  by 
Rutishauser;  that  is,  each  iterate  must  come  from  the  appropriate 
Krylov  subspace  associated  with  the  matrix  A.  Further,  we  insist  that 
each  iterate  be  optimal  over  the  Krylov  subspace  with  respect  to  a  norm 
associated  with  the  inner  product.  (The  inner  product  may  depend  upon 
A.)  Under  these  assumptions  it  can  be  shown  that  an  s-term  conjugate 
gradient  method  exists  for  every  initial  error  if,  and  only  if, 

1)  the  minimal  polynomial  of  A  is  of  degree  s 

2)  A  is  normal  with  respect  to  the  inner  product  and  A*  is 
a  polynomial  in  A  of  degree  _<  s-2  in  any  orthogonal 
basis  imposed  by  that  inner  product. 

This  implies  that  a  3-term  conjugate  gradient  iteration  exists  if, 
and  only  if.  A*  =  A  or  A  =  dl-B,  B*  =  -B,  or  the  degree  of  the  minimal 
polynomial  of  A  is  <  3. 

This  talk  will  sketch  a  proof  of  these  results  and  discuss 
implications  and  extensions. 


Session  4A.  7:30  -  7:50  p.m. ,  Monday 


A  FAST  ALGORITHM 

THAT  MAKES  MATRICES  OPTIMALLY  SPARSE 


S.  Thomas  McCormick 
Stanford  University 
Stanford,  California  94305 


Under  a  non-degeneracy  assumption  on  the  non-zero  entries  of  a 
given  sparse  matrix,  a  polynomial ly-bounded  algorithm  is  presented  that 
performs  row  operations  on  the  given  matrix  which  reduce  it  to  the 
sparsest  possble  matrix  with  the  same  row  space.  For  each  row  of  the 
matrix,  the  algorithm  performs  a  maximum  cardinality  matching  in  the 
bipartite  graph  associated  with  a  submatrix  which  is  induced  by  that 
row.  The  dual  of  the  optimal  matching  then  specifies  the  row 
operations  that  will  be  done  on  that  row.  A  variant  of  the  algorithm 
that  processes  the  matrix  in  place  is  also  described.  A  particularly 
promising  application  of  this  algorithm  is  to  reduce  linear  constraint 
matrices  as  a  way  of  accelerating  optimization. 
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NULL-SPACE  METHODS  FOR  LARGE-SCALE  QUADRATIC  PROGRAMMING 


Walter  Murray 
Stanford  University 
Stanford,  California  94305 


Null -space  methods  form  a  powerful  class  of  techniques  for  hotn 
convex  and  non-convex  quadratic  proyrammi ng ,  We  shall  consider  such 
methods  for  the  solution  of  two  classes  of  large  sparse  QP.  The  first 
class  contains  problems  for  which  t,  the  number  of  constraints  active 
at  the  solution,  is  close  to  n,  the  number  of  variables.  The  second 
class  contains  problems  for  which  t  is  small  relative  to  n.  A  dominant 
feature  of  the  storage  and  computational  overhead  required  for  a  null- 
space  method  is  the  solution  of  an  (n-t)  *  (n-t)  system  of  linear 
equations.  It  will  be  shown  how  such  systems  may  be  solved  efficiently 
when  t  i s  sma 1 1 . 


Session  12B.  11:00  -  11:20  a.m.,  Wednesday 
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A  PORTABLE  SOFTWARE  PACKAGE  FOR 
NONLINEAR  PARTIAL  DIFFERENTIAL  EQUATIONS 


J.  W.  Neuberyer 
North  Texas  State  University 
Denton,  Texas  76203 

R.  J.  Renka* 

Union  Carbide  Corporation,  Nuclear  Division 
Oak  Ridge,  Tennessee  37830 


This  paper  describes  a  storage-efficient  method  and  associated 
software  package  for  the  solution  of  a  general  class  of  second-order 
nonlinear  partial  differential  equations.  The  method  is  essentially 
an  iterative  scheme  for  the  least  squares  minimization  of  a  finite 
difference  discretization  of  the  residual.  The  advantages  of  this 
approach  are  its  ability  to  treat  a  wide  range  of  problems  (it  is 
type-independent)  and  its  low  storage  requi  rements  -  0(N)  for  N  grid 
points.  The  effectiveness  of  the  method  is  based  on  a  gradient- 
smoothing  technique  which  increases  the  rate  of  convergence  of  the 
iterative  procedure. 


Session  2B.  11:40  -  12:00  noon,  Monday 


ELEMENT  PRECONDITIONING 


Bahrain  Nour-Omid* 
Beresford  N.  Parlett 
University  of  California 
Berkeley,  California  94720 


Elliptic  boundary  value  problems  are  turned  into  large  symmetric 
systems  of  equations.  These  large  matrices  are  usually  assemDled  from 
small  ones.  It  Is  simple  to  omit  the  assembly  process  and  use  tne  code 
to  accumulate  the  product  Kv  for  any  v.  Consequently  the  conjugate 
gradient  algorithm  (CG)  can  be  used  to  solve  Ku  =  f  without  ever 
forming  K. 

It  is  well  known  however  that  CG  should  be  applied  to  precondi¬ 
tioned  systems.  In  this  paper  we  show  how  to  achieve  preconditioning 
without  forming  any  large  matrices. 

The  trade  off  between  time  and  storage  is  examined  for  the  2-3 
model  problem  and  the  analysis  of  several  realistic  structures. 


Session  12A.  11:40  -  12:00  noon,  Wednesday 


BANDED  PRECONDITIONING  FOR  SOLUTION  OF 
FINITE  ELEMENT  EQUATIONS 


Bahram  Nour-Omid* 
Horst  D.  Simon 
University  of  California 
Berkeley,  California  94720 


The  preconditioned  conjugate  gradient  algorithm  has  been 
succesful ly  applied  to  solving  symmetric  linear  systems  of  equations 
arising  from  finite  difference  and  finite  element  di scret i zat i ons  of 
a  variety  of  problems. 

In  this  paper  we  consider  a  matrix  splitting,  A  -  M  -  R,  where  M 
is  tne  dense  banded  part  of  A.  We  identify  a  certain  class  of 
symmetric  positive  definite  matrices  for  which  M  is  also  positive 
definite.  This  class  contains  matrices  arisiny  from  finite  element 
discretization  of  elliptic  boundary  value  problems. 

M  was  used  as  a  preconditioning  matrix  for  the  conjugate  gradient 
and  the  Lanczos  algorithms  to  solve  a  variety  of  practical  problems  in 
structural  engineering.  Favorable  results  were  obtained. 


Session  12A.  10:40  -  11:00  a.m.,  Wednesday 


SPARSE.  MATRICES  IN  SIMPLIFICATION  OF  SYSTEMS 


M  .  S  •  K«.  <3  » ) 

On  we  rb  ’  r.y  of  Wisconsin  -  Pawside 
Kenosna,  Wisconsin  63140 


In  recent  years  simplification  of  hi  oner  order  systems  has  been 
investigated  by  several  authors  and  the  simplified  models  are  used  for 
various  design  purposes.  The  controller  design  using  simplified  models 
for  computer  control  application  requires  that  the  simplification  algo¬ 
rithm  employed  to  derive  lower  order  models  use  minimal  computer  tine 
and  memory.  It  is  known  that  the  computat i onal  requirements  of  tne 
Routh  approximant  simplification  procedure  are  minimal. 

To  derive  lower  order  models  for  stable  systems  the  Routh  approxi¬ 
mant  method  involves  transf ormat i on  of  the  system  description  into  v- c 
canonic  form  and  then  suppression  of  less  dominant  'nodes.  The  simpli¬ 
fication  of  unstable  systems  using  Routh  approximant  procedure  involves 
a  translation  equivalent  to  the  shift  of  imaginary  axis  in  the  s  plane, 
followed  by  a  transformation  of  the  modified  system  into  y- 1  canonic 
form. 


In  both  cases  the  major  computat i ona 1  requirements  are  involved  in 
the  evaluation  and  inversion  of  the  appropriate  transf ormat ion 
matrices,  which  are  sparse  and  have  special  structure.  The  paper 
presents  the  details  of  these  sparse  matrices  and  a  procedure  to  deter¬ 
mine  their  inverses. 


Session  313.  l?:bO  -  1:10  p.m.  ,  Monday 


RECONDITIONING  STRATEGIES  TOR  THE 
CONJUGATE  GRADIENT  ALGORITHM  ON  MULTIPROCESSORS 


A.  Sumeh* 

C.  Taft 

Uni  vers i ty  of  Illinois 
Urbana,  Illinois  61801-2987 


In  this  paper  we  consider  the  conjugate  gradient  (C.G.)  algorithm, 
with  preconditioning,  for  solving  positive-definite  linear  systems  on 
multiprocessors.  We  deal  primarily  with  those  systems  that  arise  from 
the  f i ni te-rfi f ference  discretization  of  self-adjoint  elliptic  boundary 
value  proDlems  in  two  and  three  dimensions.  Assuming  that  the  multi¬ 
processor  consists  of  a  set  of  linearly  connected  processors,  we 
investigate  the  organization  of  the  C.G.  algorithm  so  as  to  achieve 


investigate  the  organization 


algorithm  so 


maximum  speedup  over  the  sequential  scheme.  Preconditioning  strategies 
that  enhance  the  convergence  of  the  C.G.  algorithm  are  also  investi¬ 
gated.  For  example,  in  the  case  of  two-dimensional  problems,  precondi¬ 
tioning  strategies  that  are  suitable  for  our  multiprocessor  include: 

(a)  block-Jacobi  splitting  in  conjunction  with  line  red-black  ordering, 

(b)  incomplete  Cholesky  f  actori  zat  i  on  (for  M-matnces)  in  conjunction 
with  point  red-black  ordering,  and  (c)  a  block  generalization  of  the 
incomplete  Cholesky  factori zation  in  conjuction  with  two  (or  more)  line 
red-black  ordering.  Furthermore,  we  explore  the  suitability  of  the 
above  preconditioned  C.G.  schemes  when  the  processors  are  not  tightly 
coupled,  i.e.,  when  the  transfer  of  one  floating-point  number  from  one 
processor  to  another  is  more  costly  than  one  arithmetic  operation. 


Session  4B.  8:10  -  8:30  p.m.  ,  Monday 


LAS02--SPARSE  SYMMETRIC  EIGENVALUE  PACKAGE 


David  S.  Scott 
University  of  Texas 
Austin,  Texas  78712 


The  LASO  package  (the  Lanczos  Algorithm  with  Selective 
Orthogonal i zati on )  became  available  through  tne  National  Energy 
Software  Center  at  the  Argonne  National  Laboratory  in  1981.  LASO  used 
EISPACK  subroutines  to  compute  a  few  eigenvalues  of  a  symmetric  band 
matrix  at  each  step  of  the  algorithm.  It  soon  became  apparent  that  the 
EISPACK  subroutines  were  not  very  efficient  when  the  band  matrix  got 
very  long  compared  to  the  width  of  the  band.  A  new  eigenvalue  solver 
was  developed  based  on  the  Rayleigh  quotient  iteration  which  overcame 
this  difficulty.  This  paper  describes  the  eignevalue  solver,  its 
i ncorporati on  into  the  package  and  the  other  modifications  which  were 
made  to  the  package  at  the  same  time. 
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IMSL  SOFTWARE  FOR  DIFFERENTIAL  EQUATIONS 
IN  ONE  SPACE  VARIABLE 


Granville  Sewell 
IMSL,  Inc. 


Houston.  Texas  770. 


Two  routines  developed  by  the  author  for  the  IMSL  Library  are 

di scussed. 

DTPTB  (Edition  8)  solves  boundary  value  problems  using  a  multiple 
shooting  technique,  with  IMSL  routine  DVERK  employed  to  solve  all 

initial  value  problems.  When  Newton's  method  is  used  to  solve  the 

resulting  system  of  simultaneous  equations,  a  Jacobian  with  a  "periodic 
band"  structure  arises. 

DPDLS  (Edition  9)  solves  a  partial  differential  equation  system  of 
the  form  ut  =  f ( x ,t ,u ,ux ,uxx)  using  a  collocation  method  and  the  method 

of  lines.  Cubic  Hermite  basis  functions  are  used  to  integrate  the 

method  of  lines  ODE  systems  y'  =  A^fft.y),  where  although  both  A  and 


ines  ODE  systems 


A  f(t,y),  where  although  both 


■j  -j  ,  \  » 

the  Jacobian  of  f  are  banded,  A'1  is  full.  A  very  simple  modification 
of  DGEAR,  which  could  be  applied  to  other  GEAR  codes  designed  only  for 
the  usual  ODE  problem  y'  =  f(t,y),  was  necessary  to  handle  this  system 
without  the  need  to  deal  with  any  full  matrices. 


Session  2B.  11:20  -  11:40  a.m. ,  Monday 


THE  LANCZOS  ALGORITHM  WITH  PARTIAL  REORTHOGONAL I ZAT I  ON 
FOR  THE  SOLUTION  OF  NONSYMMETRIC  LINEAR  SYSTEMS 


Horst  D.  Simon 
SONY  at  Stony  Brook 
Stony  Brook,  New  York  11794 


The  Lanczos  algorithm  with  partial  reorthogonal i zati on  has 
recently  been  applied  by  the  author  to  the  solution  of  large  sparse 
symmetric  linear  systems  of  equations.  Here  this  work  is  extended  for 
the  nonsymmetric  case  and  linear  least  squares  problems  Bx  =  f.  This 
is  done  by  applying  the  symmetric  algorithm  to  the  problem 


as  Paige  [TOMS  8,  Vol.  1,  1982]  already  showed,  several  simplifications 
result  because  of  the  special  structure  of  the  matrix  and  the  right 
hand  side.  By  using  partial  reorthogona 1 i zat i on  a  new  algorithm  is 
obtained  which  minimizes  the  number  of  matrix-vector  multiplications. 
We  present  some  numerical  examples,  which  show  that  the  algorithm  works 
well. 
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SOME  OPERATIONS  ON  SPARSE  MATRICES 


Daniel  B.  Szyld* 

Oleg  Vi shnepol sky 
Institute  for  Economic  Analyst 
New  York  University 
New  York,  New  Yorx  1 000 3 


In  many  appl i cati ons ,  the  matrices  are  so  large  that  even  with 


sparse  matrix  storage  scheme,  not  more  than  one  or  sometimes  two  of 
them  will  fit  in  core.  We  developed  an  algorithm  for  permuted 
transposition  when  only  the  original  matrix  is  kept  in  core  and  the 
resulting  one  is  written  on  secondary  storage.  This  algorithm  is  based 
on  the  distribution  of  all  non-zero  elements  among  the  column  queues, 
using  a  link  vector. 


A  second  algorithm,  to  perform  the  multiplication  of  two  sparse 
matrices,  overwrites  the  portion  of  the  factors  that  has  already  been 
used,  and  if  this  amount  of  storage  is  not  sufficient,  it  writes  the 
product  out  of  core.  These  two  algorithms  extend  those  presented  by 
F.  Gustavson  [TOMS ,  4,  1973,  pp.  250-263].  We  present  some 
comparisons.  . . 

To  obtain  rows  and/or  columns  permutations  of  a  sparse  matrix,  one 
can  execute  the  permuted  transposi t i on  algorithm  twice.  We  offer  an 
alternative  method  that  is  faster  if  the  resulting  matrix  is  used  with 
its  indices  unordered,  and  we  explore  different  sorting  techniques  for 
ordering  the  indices.  The  same  algorithm  is  used  to  obtain  a  submatrix 
(any  subset  :f  rows  and  columns)  of  a  sparse  matrix. 
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This  talk 

will  describe  a  technique  that  permits  t.ne  representa- 

-  '»■ 

tion  of  the  inverse  of  a  matrix  A  with  only 

one  triangular  array.  More 

precisely,  let 

LA  =  U  with  L  lower  triangular  and  U  upper  triangular 

arrays.  Then, 

to  compute  A" *b  and  b  A" 1  in 

0(N2)  operations,  we  need 

the  arrays: 

_'V- 

b 

(1)  L,  U 

(Gauss) ; 

(2)  L-1, 

U  (Crout); 

(3)  L,  A 

( Imp  1 icit  Gauss) . 

i 

The  array 

L  can  be  obtained,  in  0 ( (V 3 ' 

operations,  without  having 

to  compute  all  the  elements  of  U. 

We  will  discuss  some  of  the  sparsity  implications  of  being  able  to 
concentrate  on  obtaining  a  sparse  array  L,  not  caring  about  the  number 
of  non-zeros  in  the  array  U. 


,v 
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